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Abstract. 

We develop a path-integral formalism to study the formation of large-scale structures in the universe. Starting from 
the equations of motion of hydrodynamics (single-stream approximation) we derive the action which describes the 
statistical properties of the density and velocity fields for Gaussian initial conditions. Then, we present large- N 
expansions (associated with a generalization to N fields or with a semi-classical expansion) of the path-integral 
defined by this action. This provides a systematic expansion for two-point functions such as the response function 
and the usual two-point correlation. We present the results of two such expansions (and related variants) at one- 
loop order for a SCDM and a ACDM cosmology. We find that the response function exhibits fast oscillations in 
the non-linear regime with an amplitude which either follows the linear prediction (for the direct steepest-descent 
scheme) or decays (for the 2PI effective action scheme). On the other hand, the correlation function agrees with the 
standard one-loop result in the quasi-linear regime and remains well-behaved in the highly non-linear regime. This 
suggests that these large- A expansions could provide a good framework to study the dynamics of gravitational 
clustering in the non-linear regime. Moreover, the use of various expansion schemes allows one to estimate their 
range of validity without the need of A^— body simulations and could provide a better accuracy in the weakly 
non- linear regime. 
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1. Introduction 

The large-scale structures we observe in the present uni- 
\ verse (such as galaxies and clusters of galaxies) have 
formed thanks to gravitational instability (Pebbles 1980) 
which amplified the small density perturbations created in 
the early universe (e.g. through quantum fluctuations gen- 
erated during an inflationary phase, Liddle & Lyth 1993). 
Moreover, the power increases at small scales as in the 
CDM model (Peebles 1982) which leads to a hierarchical 
scenario where small scales become non-linear first. Then, 
small structures merge in the course of time to build in- 
creasingly large objects. Therefore, at large scales or at 
early times one can use a perturbative approach to de- 
scribe the density field (e.g. Bernardeau et al. 2002 for 
a review). This is of great practical interest for observa- 
tions such as weak-lensing surveys (e.g. Bernardeau et al. 
1997) and CMB studies (e.g. Seljak & Zaldarriaga 1996). 
At smaller scales one uses N-body numerical simulations 
to obtain fits to the density power-spectrum (e.g. Smith 
et al. 2003) which can be compared with observations. 

In the weakly non-linear regime one may apply the 
standard perturbative expansion over powers of the ini- 
tial density fluctuations (e.g., Fry 1984; Goroff et al. 1986). 
This procedure is used within the hydrodynamical frame- 
work where the system is described by density and ve- 



locity fields rather than by the phase-space distribution 
/(x, v, t), that is one starts from the Euler equation of mo- 
tion rather than the Vlasov equation. This is quite success- 
ful at tree-order where one can obtain the leading-order 
contribution to all p— point correlations and reconstruct 
the density probability distribution (Bernardeau 1994), 
assuming Gaussian initial fluctuations as in usual infla- 
tionary scenarios (Liddle & Lyth 1993). However, several 
problems show up when one tries to compute higher-order 
corrections. First, ultraviolet divergences appear for linear 
power-spectra Pi(fc) oc k n with n > — 1 (Scoccimarro & 
Frieman 1996b) which would break the self-similarity seen 
in numerical simulations. This may also be interpreted as a 
consequence of the breakdown of the hydrodynamical de- 
scription beyond shell-crossing (Valageas 2002). However, 
for n < — 1 the one-loop correction yields a good agree- 
ment with numerical simulation (Scoccimarro & Frieman 
1996b). Secondly, being an expansion over powers of Pb(k) 
this standard perturbative approach yields a series of 
terms which grow increasingly fast into the non-linear 
regime at higher orders and cannot relax at small scales. 
This makes the series badly behaved and it seems difficult 
to go deeper into the non-linear regime in this manner. 

On the other hand, a good description of the weakly 
non-linear regime becomes of great practical interest as 
cosmological probes now aim to constrain cosmological pa- 
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rameters with an accuracy of the order of 1% to keep pace 
with CMB observations. In particular, weak-lensing sur- 
veys mainly probe this transition regime to non-linearity 
and a good accuracy for the matter power-spectrum is 
required to make the most of observations and constrain 
the dark energy equation of state (Munshi et al. 2007). 
Another probe of the expansion history of the Universe 
is provided by the measure of the baryon acoustic oscil- 
lations which also requires a good description of weakly 
non-linear effects (Seo & Eisenstein 2007; Koehler et al. 
2007). 

As a consequence, it is worth investigating other ap- 
proaches which may partly cure some of these problems 
and provide a better description of the non-linear regime. 
For instance, McDonald (2006) proposed a renormaliza- 
tion group method to improve the standard perturba- 
tive approach. On the other hand, Crocce & Scoccimarro 
(2006a, b) presented a diagrammatic technique to organize 
the various terms which arise in the standard perturbative 
expansion and to perform some infinite partial resumma- 
tions. Besides, using the asymptotic behavior of the ver- 
tices of the theory they managed to complete the resum- 
mation in the small-scale limit for the response function 
R (also called the propagator). 

In this article, we apply to the hydrodynamical frame- 
work the method presented in Valageas (2004) for the 
Vlasov equation. This approach introduces a path-integral 
formalism to derive the statistical properties of the system 
from its action S. Then, one applies a "large-A^ expan- 
sion" (which can be seen as similar to the semi-classical 
expansion over powers of h or can be derived from a gener- 
alization to N fields) to compute the quantities of interest 
such as the two-point correlation. This offers a different 
perspective to address the non-linear regime of gravita- 
tional clustering which complements other approaches and 
may also serve as a basis for other approximation schemes 
than the ones described in this article. For instance, previ- 
ous approaches can be recovered by expanding the path- 
integral obtained in this paper in different manners (e.g. 
the usual perturbative expansion corresponds to the ex- 
pansion over the non-Gaussian part of the action S) . The 
two large- N expansions discussed in this article are two 
other means of performing systematic expansions to com- 
pute the required correlation functions. This also amounts 
to reorganize the standard perturbative expansion by per- 
forming some partial infinite resummations. In this fash- 
ion, one can hope to obtain new expansion schemes which 
are better suited to describe the non-linear regime of grav- 
itational clustering. 

First, we recall in sect. [2] the equations of motion ob- 
tained within the hydrodynamical framework. Next, we 
derive the path-integral formulation of the problem in 
sect. [3] and we describe two possible large- N expansions 
in sect. 2] Finally, we present our numerical results in 
sects. [5][7l focusing on a critical-density universe. Since 
our formalism applies equally well to any cosmology we de- 
scribe our results for a ACDM universe in sect. [8] Finally, 
we conclude in sect. [9] 



2. Equations of motion 

2.1. Fluid approach 

At scales much larger than the Jeans length both the cold 
dark matter and the baryons can be described as a pres- 
sureless dust. Moreover, before orbit crossing one can use 
a hydrodynamical approach (in the continuum limit where 
the mass m of particles goes to zero) . Then, at scales much 
smaller than the horizon where the Newtonian approxima- 
tion is valid, the equations of motion read (Peebles 1980): 

^ + V.[(l + £)v]=0, 



dv 

— + Hv + (v.V)v = -V(/>, 

OT 



= -n m n 2 s, 



(i) 

(2) 
(3) 



where r = J dt/a is the conformal time (and a the scale 
factor), 7i = dlna/dr the conformal expansion rate and 
f2 m the matter density cosmological parameter. Here S 
is the matter density contrast and v the peculiar veloc- 
ity. Since the vorticity field decays within linear theory 
(Peebles 1980) we take the velocity to be a potential field 
so that v is fully specified by its divergence 9 = V.v. Then, 
in Fourier space eqs.fl])-([3]) read (Goroff et al. 1986): 

<9<5(k,T 



Or 



0(k, t) = / dkidk 2 fo(ki + k 2 - k) 
xa(k 1 ,k 2 )0(k 1 ,r)5(k 2 ,r), 



(4) 



&l)+^(k,r) + ^Q ni W 2 J(k,r) = 
or 2 

-J dkxdk 2 <J D (ki + k 2 - k)/J(ki,k a )0(ki,T)0(k a ,T), (5) 

where 8d is the Dirac distribution, the coupling functions 
a and /3 are given by: 

(ki+k 2 ).ki |ki+k 2 | 2 (ki.k 2 ) 
a(k 1; k 2 ) = -j ,/3(ki,k 2 ) = — ^ , (6) 

and we defined the Fourier transforms as: 
dx 



5(k) 



-2k. x 



(2tt) s 



a(x) 



(7) 



As in Crocce & Scoccimarro (2006a, b) let us define the 
two-component vector ip as: 



ib(k n)- ( _ ( SQt,tj) \ 

~~ ' MKv) J ~ {-e(k,rj)/Hf J 



(8) 



where we introduced the time coordinate 77 defined from 
the linear growing rate D + of the density contrast: 



, D+(t) , dlnD+ 1 dln£>4 

v = 1 ^— — , / 



(9) 



D + q din a Ti dr 

with D + o is the value of D + today at redshift z — and: 

-n m n 2 D+. (10) 



d 2 D+ . „.d£>+ 3. 



dr 2 



n- 



dr 
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Then, the equations of motion (@|-([5]) can be written as: 

0(x,x').ip(x') = K a (x;xi,X2).ip(xi)ip(x2), (11) 

where we introduced the coordinate x = (k, 77, i) where 
i — 1, 2 is the discrete index of the two-component vectors. 
In eq. (jTTJ) and in the following we use the convention that 
repeated coordinates are integrated over as: 

0(x, x')4{x') = J dk'dr?' <V( k , 77; k', r/)Vi'(k', »/)-(12) 

i'=l 

The matrix O reads: 

0(x,x')=( _I m b_ , ~sL _ 1 )s D (k-k , )S D (r 1 -r 1 ')(13) 

whereas the symmetric vertex K s (x; x\, x 2 ) = 
K s (x]X2,x\) writes: 

K s (x;xi,x 2 ) = S D (ki + k 2 - k)5 D (r)i - r))5 D (r}2 - 77) 

x 7 ? ilii2 (k 1 ,k 2 ) (14) 



Thus, for a CDM cosmology the linear power A L0 (k) 
grows as k 4 at small k and as lnfc at high k. Then, the 
linear two-point correlation function Gl{x\ 1 x 2 ) reads: 



Gl{xi,x 2 ) = (iPl(xi)^l(x2)) 



= e" 1+ " 2 F L0 (fc 1 )fe(k 1 +k 2 ) 



1 1 
1 1 



(21) 



As in Valageas (2004) it is convenient to introduce the 
response function R(x\, x 2 ) (which is also called the prop- 
agator in Crocce & Scoccimarro (2006a, b)) defined by the 
functional derivative: 



R(xi,x 2 ) 



,5ip(xi) 



>C=o, 



(22) 



with: 

7l:l, 2 ( k l, k 2) 



a(k 2 ,ki) a(ki,k 2 ) , , 

> 7l;2.l( k l> k 2) = o ,(15) 



2 ' " 2 

7 f. 2>2 (k 1) k 2 )=/3(k 1 ,k 2 ), (16) 

and zero otherwise (Crocce & Scoccimarro 2006a). 
2.2. Linear regime 

At large scales or at early times where the density and ve- 
locity fluctuations are small one can linearize the equation 
of motion (Ti"Tj) which yields O.ipL = 0. Then, the linear 
growing mode is merely: 



M*) =e"<J L0 (k) 



(17) 



where <Sz,o(k) is the linear density contrast today at red- 
shift z = 0. We shall define the initial conditions by the 
linear growing mode (|17[) (which is consistent with the 
fact that we only kept the curl-free part of the peculiar 
velocity in eqs.(|l])-([5])). Moreover, we assume Gaussian ho- 
mogeneous and isotropic initial conditions defined by the 
linear power-spectrum PLo(k): 



(<M k i)<M k 2)) = ^o(fci)M k i 



(18) 



It is convenient to define the power per logarithmic 
wavenumber A 2 (fc) by: 



A 2 (k) = 47rfc 3 P(fc), 
whence: 

. , . , f°° dk 2 sinfc|xi-x 2 | 

(6( X1 )5(*2)) = J o T A{k) fc|xi _ x2| 



(19) 



(20) 



S((X2) ' 

where ((x) is a "noise" added to the r.h.s. in cq. lfTTj) . Thus, 
R(xi, x 2 ) measures the linear response of the system to an 
external source of noise. Because of causality it contains 
an Heaviside factor #(771 — 772) since the field ip(xi) can 
only depend on the values of the "noise" at earlier times 
V2 5: *7i- Moreover, it satisfies the initial condition: 



77^ : R(xi,x 2 ) -> S D (ki - k 2 )5 ilt 



(23) 



In the linear regime, the response Rl may be obtained 
by first calculating the linear solution tpL{x\C,) of the lin- 
earized equations of motion 0.iPl{x\C) = C an d next tak- 
ing the functional derivative (|22"|) . Alternatively, one can 
directly obtain the matrix Rl from the initial condition 
(|23|) and the linear dynamics O.Rl — (as implied by the 
definition (f2"2")) and O.ipL = 0). For the Einstein-de-Sitter 
cosmology (O ro = 1, Oa = 0) where the factors f2 m // 2 are 
constant and equal to unity, this can be easily solved by 
looking for eigenmodes or performing a Laplace transform 
as in Crocce & Scoccimarro (2006a). This yields: 



Rl(xi,x 2 ) = 0(771 - 77 2 )<5_D(ki - k 2 ) 
3 2 



2 -2 
-3 3 



(24) 



For other cosmological parameters one must solve numer- 
ically the linearized equation of motion for Rl . Note that 
in all cases it exhibits the factor 6(7)1— ^2)^D(ki— k 2 ) times 
a term which is independent of wavenumbers ki, k 2 . 

3. Path-integral formalism 

As in Valageas (2004) we can apply a path-integral ap- 
proach to the hydrodynamical system pip since we are 
only interested in the statistical properties of the density 
and velocity fields (we do not look for peculiar solutions 
of the equations of motion). Let us briefly recall how this 
can be done (also Martin et al. 1973; Phythian 1977) . In 
order to include explicitly the initial conditions we rewrite 
cq. ([TTj) as: 



with ip = for 77 < rji and: 



(25) 



(26) 
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Thus, the source fii merely provides the initial conditions 
at time rji, obtained from the linear growing mode (|17[) . 
We shall eventually take the limit rji — ► — oo. Next, we 
define the generating functional Z[j] by: 



m = 



= / [dMi] e 



(27) 



where we took the average over the Gaussian initial con- 
ditions: 



Thus, the statistical properties of the system (JTTJ) are de- 
scribed by the action S[ip, A] defined by: 

S[ip, A] = X.(0.ip - K,.iM>) - ^A.Aj.A. (36) 

Moreover, we can note that adding a "noise" £ to the r.h.s. 
of eq. amounts to change fa — ► fj,i + £ which translates 
into 5 — > 5 — A.C- Therefore, functional derivatives with 
respect to C are equivalent to insertions of the ghost field 
A. In particular, we have: 



(/it) = 0, (fJ, i (x 1 )fX i (X2)) = A i (x 1 ,x 2 ). 



(28) 



All statistical properties of the field ip may be obtained 
from Z[j}. It is convenient to write eq. (|27|) as: 



R{xx,X2)={ip(x 1 )X{x 2 )), (A>=0, (AA>=0. 



4. Large-A" expansions 



m = 



xe 



[dM i ][dV'] |detAf| fe(^ 



04 + KM) 



(37) 



The path-integral (f3"5| can be computed by expanding 
over powers of its non-Gaussian part (i.e. over powers of 
K s ). This actually yields the usual perturbative expan- 
sion over powers of the linear power-spectrum Pl (see 
also Valageas (2001, 2004) for the case of the Vlasov equa- 
tion of motion). On the other hand, the path-integral (|35[) 
may also be studied through a large- TV expansion as in 
Valageas (2004). Thus, one considers the generating func- 
detM= det(0-2K s ,ip) = det(d/d V ) det(l+O-2K a .ip)(30) tional Z N [j,h] defined by: 



(29) 



where the Jacobian |detM| is defined by the functional 
derivative M = Sfii/Sip. We can compute the latter by 
factorizing the operator d/drj: 





3^(ry') 



-1 
M^O/) 



with: 

d(x,x')=S D (\i-k')9(ri-ri') 

2/ 2 VI ) 2/ 

(note that we separated the identity part out of O) and: 

K s (x; xi,x 2 ) = S D (ki + k 2 - k)9(r) - 771)5.0(772 ~ m) 

x 7 |. ili , 2 (k 1 ,k 2 ). (32) 

Then, we have: 

det M = det(d/d V ) e Trln(l+0-2K s .^) 



Z N \j,h}= /[d#dA] 



N\j.^+h.\-S[1>,\]] 



(38) 



and: 

Trln(l + 



2K s .i>) 



(-1) 



9 -l 



1=1 



-Tr(<5 - 2K s .iP) q 



= Tr(G - 2K S 4), 



where we used the fact that the trace of powers q > 2 
vanishes because of the Heaviside factors 9(r) — 77') (Zinn- 
Justin 1989, Valageas 2004). Next, the trace TrO is an ir- 
relevant constant whereas Tr(2K s .ip) = (we regularize as 
usual the gravitational interaction at large distances to re- 
move ambiguities, e.g. perform angular integrals first, see 
also Valageas 2004). Therefore, the Jacobian in the func- 
tional integral (f2"9"| is a constant which can be absorbed 
into the normalization. Then, introducing an imaginary 
ghost field A to express the Dirac as an exponential and 
performing the Gaussian integration over [ii we obtain: 



Z\j] = /[d/i^pA] e 



j.ip+x.(m-o.ip+K B .ipiii)-ifj,i.Ar 1 ,n i 



1 j (31) 

and one looks for an expansion over powers of 1/A, tak- 
ing eventually A — 1 into the results. As discussed in 
Valageas (2004) the large-A expansions may also be de- 
rived from a generalization of the dynamics to A fields 
V> (q) . This yields the same results once we take care of the 
long-wavelength divergences which constrain which sub- 
sets of diagrams need to be gathered. 

(33) 4.1. Steepest-descent method 

A first approach to handle the large-A limit of eq. (|38|) is to 
use a steepest-descent method (also called a semi-classical 
or loopwise expansion in the case of usual Quantum field 
theory with Ti — 1 /A) . This yields for auxiliary correlation 
and response functions Go and Rq the equations (Valageas 
2004): 

O(x,z).G (z,y) = (39) 
0(x,z).Ro(z,y) = S D (x~y) (40) 
Ro(x,z).0(z,y) = 5 D (x-y) (41) 

whereas the actual correlation and response functions 
obey: 

0(x,z).G(z,y) = X(x,z).G(z,y) + IL(x,z).R T (z,y) (42) 
0(x, z).R(z, y) = 5 D (x - y) + z).R(z, y) (43) 
R(x, z).0(z, y) = S D (x~y)+ R(x, z).E(z, y) (44) 

where the self-energy terms £ and II are given at one-loop 
order by: 



(34) 



= / [dr))][d\] e J-*+A.(-0.*+X..iW)+iA.A,.A i 



Y,{x, y)=AK s {x; xi,x 2 )K s (z; y, z 2 )R (xi, z)G (x 2 ,z 2 ) (45) 
( 35 ) IL(x, y)=2K s (x; Xl ,x 2 )K s (y; y 1 ,y 2 )G (x 1 ,y 1 )G (x 2 ,y 2 lA6) 
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Note that eas. (j59")) - (l4"4l are exact and that the expansion 
over powers of 1 /N only enters the expression of the self- 
energy (14"S"1) - (|4"6"1) . Here we only kept the lowest-order terms 
(see Valageas (2004) for the next-order terms). We also 
took the limit rji ~ * — oo so that terms involving A, van- 
ish. From sect. !2.2l and eas.(|39 p -([4i p we find that the aux- 
iliary matrices Go and Rq are actually equal to their linear 
counterparts: 

G = G L , R = Rl- (47) 

4.2. 2PI effective action method 

As described in Valageas (2004) a second approach is to 
first introduce the double Legendre transform T [i/j, G] of 
the functional W = h\Z (with respect to both the field 
tp and its two-point correlation G) and next to apply the 
1/N expansion to T. This "2PI effective action" method 
yields the same equations (|4"2")) - (l4"4"|) and the self-energy 
shows the same structure as (|45[) -(|46 |) where Go and i?o 
are replaced by G and R. Thus, the direct steepest-descent 
method yields a series of linear equations which can be 
solved directly whereas the 2PI effective action method 
gives a system of non-linear equations (through the de- 
pendence on G and R of E and n) which must usually 
be solved numerically by an iterative scheme. However, 
thanks to the Heaviside factors appearing in the response 
R and the self-energy E these equations can be solved 
directly by integrating forward over time 771. 

4.3. General properties 

As discussed in Valageas (2004) both the steepest-descent 
and 2PI effective action methods agree with the standard 
perturbative analysis over powers of Pho{k) up to the or- 
der used for the self-energy (e.g. up to order P£ if we 
only consider the one-loop terms (|45|) - (|46]) ). As compared 
with the standard perturbative approach, the two schemes 
described above also include two different infinite partial 
resummations, as can be seen from eqs.([4"2" j) - (|4"4")) which 
clearly generate terms at all orders over P^o for G and R 
even if E and II are only linear or quadratic over Plo- 

We can note that the equations we obtain for the hy- 
drodynamical system ([l}-([3]) are simpler than for the col- 
lisionless system studied in Valageas (2004) which is de- 
scribed by the Vlasov-Poisson system. Indeed, here the 
mean (ip) vanishes at all orders. This can be explicitly 
checked from eqs.ld])-© and in the derivation of eqs. (j39")) - 
(|46|) . This is not the case for the Vlasov dynamics where 
eq.([47|) does not hold. 

Using a diagrammatic technique Crocce & 
Scoccimarro (2006a, b) derived equations l|42 p -(|46 jl 
in an integral form (i.e. without the differential operator 
O in the l.h.s.) by first integrating the equation of motion 
(fTTj) . Of course, our path-integral approach can also be 
applied to the integral form of eq. fTTj) . This amounts to 
write cq. (|I 1| as ip = ipL + Ks^^P where tpL is the linear 
growing mode and K s the integral form of the vertex 



K s . Then, the procedure presented in sect. [3] can be 
applied to this integral form of the equation of motion, as 
described for instance in Valageas (2001). In this article 
we preferred to keep eq. (jTTJ) in its differential form so that 
the r.h.s. does not contain an integral over time. Then the 
self-energy terms £(771,172) and n(?7i , 772) only depend on 
the response and correlation at the same times (771,772), 
see eas. p5t - p6|) . By contrast, the integrated vertex K s 
would lead to self-energies which depend on the values 
of the response and correlation at all past times which 
would entail less efficient numerical computations. 

4.4. Integral expression for the two-point correlation 

Thanks to statistical homogeneity and isotropy the matri- 
ces Go , G and II are symmetric and of the form: 

G(xi,x 2 ) = S D (ki +k 2 )Gi li i 2 (fci;77i,77 2 ), (48) 

with: 

G iui2 (k; 771,773) = G i2lil (fc; 772, 77i), (49) 

whereas the matrices Rq , R and E are of the form: 

R(xi,x 2 ) = S D (ki - k 2 )0(?7i - 772)-Ri 1 ,i 2 (fei;77i,77 2 ). (50) 

Besides, from cqs. (l4"2"1) - (|4"4"]) we see that the correlation G 
can also be written in terms of the response R as: 

G{xi,x 2 ) = R x Go{rn) x R T + R.Il.R T (51) 

where the second product is defined as in (fT2|l whereas the 
first product does not contain any integration over time: 

R x G o (?70 x R T = 8 D (ki +k 2 ) Ri 1 ,i> 1 (ki;T)i,r)i) 

xGo;i' 1 ,i' 2 {ki;r)i,'ni) R i2,i' 2 ( k i;m, Vi), (52) 

and we let rji — > —00. The physical meaning of eq. (|5ip is 
clear. The first term in the r.h.s. means that the fluctua- 
tions at the initial time 77, are merely transported forward 
in time through the matrix R. This is the only non-zero 
term in the linear regime (with R = Rl and Go = Gl). 
The effect of the non-linear dynamics is to modify the 
transport matrix R and to add a second term to the r.h.s. 
of eq. ([5Tj) . The latter has the meaning of a source term 
which produces fluctuations with two-point correlation 
11(77^,772) at the times (771,772) which are next transported 
forward to later times (771,772) by the matrices R{rji,rji) 
and i? T (?7 2 , 772). Thus, the "dot product" means that we 
need to integrate over all fluctuations produced at all ear- 
lier times rj^ < rji and rj' 2 < 772. This interpretation con- 
siders II as an external source for an "effective" linear dy- 
namics. For the system which we study here the "source" 
II is not external but corresponds to the non-linear inter- 
action of the field tp, as seen from the explicit expression 
(|46l) . Moreover, the linear operator which would define the 
transport with time also depends on G and R through E, 
see eqs.(|i5|) and (|34]l . 
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Alternatively, Crocce & Scoccimarro (2006a) pointed 
out that the expression (|51[) is somewhat similar to the 
result obtained within a phenomenological halo model 
(Seljak 2000). Indeed, the latter model splits the non- 
linear power-spectrum into two terms, one that dominates 
in the linear regime (2— halo term) and the other that 
dominates in the highly non-linear regime (1— halo term). 
They could be identified with the first and second terms 
of eq. (|5"Tj) . However, it is not clear whether the analogy 
can be pursued much further. In particular, eq. (|51|) is ex- 
plicitly dynamical (i.e. it explicitly involves integrations 
over past events) whereas the halo model gives a static 
expression (it writes the two-point correlation as integrals 
over the current halo distribution, which is the basic time- 
dependent quantity). 

Note that within the 2PI effective action method, even 
if R would be treated independently of G the expres- 
sion (j5"Tj) does not give a quadratic dependence of G on 
R since LT depends quadratically on G. Besides, solving 
eq. (j5"Tj) perturbatively over Go at fixed R actually gener- 
ates terms at all orders q over powers Gq. By contrast, in 
the steepest-descent approach where n only depends on 
Go the expression (f5Tj) is fully explicit and quadratic over 
both Go(= Gl) and R if the latter are treated indepen- 
dently. This is a signature of the additional resummations 
involved in the 2PI effective action method as compared 
with the direct steepest-descent approach. 

An advantage of eq. (f5Tj) is that once II and R are 
known, we obtain an explicit expression for G. Hence we 
do not need to solve eq. (j42|) (by solving for R we have 
already performed the "inversion" of O — E). A second 
advantage of eq. (15ip is that it provides an expression for 
G which is clearly symmetric. Moreover, since we start 
from a two-point correlation Go(= Gl) which is positive, 
as shown by ea. (|2"Tj) , we see from eq. (|46|) and eq. (f5"Tj) that 
both n and G are positive. This also holds for the 2PI 
effective action approach which may be obtained by iter- 
ating the system (j42 |) - (|46)) (substituting G and R into the 
self-energy). 

4.5. Integration over wavenumber for the self-energy 

Using the symmetry (|49|) for n and the form (|50)) for E we 
see that H(x\, x 2 ) and E(xi, X2) only need to be computed 
at times 771 > r\i- Besides, thanks to the various Dirac 
factors the expressions ([4"5"j) -([4l) |) only involve an integral 
over one wavenumber k'. Thus, using the expression (|14p 
of the kernel K s we can write from eqs. (|45l) - (|46|) the self- 
energy as: 

E ni2 (fc; m,m) = 4 7 ? A)j2 (k',k - k') 7 f i;W2 (k, k' - k) 
xR h i 1 (k';r] 1 ,T]2)Gj 2 i 2 {\k- k'|; 771, 772) (53) 



and: 



U ili2 (k; Vl ,V2) = 2 7 « y . ij2 (k',k-k07^ ;il , i2 (-k / ,-k + k') 
xG nh (k'; rn,ri 2 )G j2 i 2 (\k - k'|; 771, 772), (54) 



where we integrate over k' and we sum over ji,j 2 , h,h- In 
order to keep the symmetry k' <-> k — k' in the numerical 
integration we use symmetric elliptic coordinates: 



k'=* 

2 



with q = — 
H 2 



sinh £ sin /_t cos < 
sinh C sin /i sin q 
cosh £ cos ji 



(55) 



where k is taken along the third axis and £ > 0, < fj, < 

7T, < (j>< 2tt. 

At equal times 771 = r\ 2 we can obtain from eq. (|53p : 



E(fc; 77,77) 



47T , 



dk'G 2 2(k';v,v) x 



1 
1 



(56) 



The explicit expression (|56|) shows that within the direct 
steepest-descent method we need a linear power-spectrum 
PLo(k) with a logarithmic slope n = d In Plo/A In k such 
that n > — 1 at large scales and n < — 1 at small scales 
to obtain finite results. These constraints are the same 
as those which appear in the standard perturbative ap- 
proach. The divergence at long wavelengths eventually 
cancels out if we only consider the equal-time correlations 
of the density field, as seen in Jain & Bertschinger (1996) 
(also Vishniac 1983). As pointed out by these authors, this 
can be understood from the fact that contributions to the 
velocity field from long-wave modes correspond to an al- 
most uniform translation of the fluid and therefore should 
not affect the growth of density structures on small scales. 
However, these infrared contributions can no longer be ig- 
nored when one studies different-times correlations or the 
velocity field itself. The same problem appears when one 
considers the collisionless dynamics as described by the 
Vlasov equation (Valageas 2004). The small-scale diver- 
gences are more harmful as they do not cancel in such a 
fashion. In fact, they may be understood as the signature 
of anomalous terms (which do not scale as integer pow- 
ers of Plo) which cannot be recovered by the standard 
perturbative scheme as discussed in Valageas (2002). It is 
probably necessary to go beyond the hydrodynamical ap- 
proximation |l|-(I3]) to cure this small-scale problem (e.g. 
use the Vlasov equation which does not break down at 
shell-crossing). On the other hand, note that within the 
2PI effective action scheme these constraints apply to the 
non-linear power G 22 (k) which may have different a differ- 
ent asymptotic behavior than the linear power-spectrum 
at small scales. Therefore, the range of initial conditions 
where the integral in eq. (|56[) converges at high k may be 
different within this scheme. 

5. Direct steepest-descent method 

We present in the following sections our numerical results 
which we compare to fits obtained from direct numeri- 
cal simulations and to the standard perturbative anal- 
ysis. We consider an Einstein-de-Sitter cosmology with 
il m = 1, — 0, a shape parameter T = 0.5 and a 
normalization eg = 0.51 for the linear power-spectrum 
Pio(fc) and a Hubble constant H = 50 km.s~ 1 .Mpc~ 1 
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(i.e. the reduced Hubble parameter is h = 0.5) as in Smith 
et al.(2003) so as to compare our results with N-body sim- 
ulations. Since we mainly wish to investigate in this arti- 
cle the properties of the large- N expansions described in 
sect. 3] we first focus on the Einstein-de-Sitter cosmology 
where both Gl and Rl have simple analytical forms (|2"Tj) - 
(|24|) . We shall discuss the case of a ACDM cosmology in 
sect. [8] below. 

We first study in this section the results obtained from 
the direct steepest-descent method of sect. 14.11 which in- 
volves the auxiliary two-point functions Go and Rq. Since 
the latter are equal to their linear counterparts (|47|) the 
correlation G is given by eq. (f2"Tj) whereas i?o is obtained 
by solving eq. (|40p forward in time, in agreement with the 
discussion below eg. ((23)) . For the Einstein-de-Sitter cos- 
mology which we consider here Rq is actually given by 
the explicit expression ([24|) . 

5.1. Self-energy 



1 oo 



10 



0.1 



0.01 



o.o : 




0. 1 1 

k [h Mpc- 1 ] 



Fig. 1. The self-energy terms Ej^ (fc) of eq.(|60|) as a func- 
tion of wavenumber k. We display the four components 
^o-n (solid line), E^ 2 i (dot-dashed line), Sj 12 (dotted 
line) and E J 22 (dashed line) . They are all negative except 
for Ej n which is positive at k < 0.2h Mpc -1 and nega- 
tive at higher k. All terms are of the same magnitude and 
grow roughly as k 2 in agreement with eq. (|62l) . 



From Go and i?o we need to compute the self-energy 
from eqs. (|45l) - (l46|) . In fact, for the Einstein-de-Sitter cos- 
mology the time-dependence of the self-energy can be fac- 
torized as follows using eqs. ([2"Tj) . ([24]) . First, the auxiliary 
response Rq can be written from eq. ((24|) as: 



R (xi,x 2 ) = 0(77i - T)2)Sd0^i -k 2 ) 



e r ' 1 - T > 2 R+ + e-^ r > 1 - r ' 2 ^ 2 R- 



(57) 



with: 



(3/5 2/5 1 i,ml '''' 



2/5 -2/5 
-3/5 3/5 



(58) 



where we separated the contributions from the growing 
mode Rq and the decaying mode Rq . Then, from eq. ([53")) 
the self-energy term £0 which is linear over i?o also splits 
into two terms as: 

£0(2:1, #2) = (r/x - r/ 2 )S D (ki -k 2 ) 

x e 2 ^£+(fc 1 ) + e-^/ 2+5f ' 2 / 2 I] (fci)] ,(59) 

with: 

^Jk) = 4 E / d <17^ lj2 (k-q,q) 



jihhl 



X 



^ ;ia , ia (k ) -q)P M (?)^ ji(l 



(60) 



where we used eq. (|2"Tj) . Besides, from the expressions ([T5 
(fT6| of the vertices 7 s we obtain: 



^0 (k) — ( _v + 2 v 

^0:21 ^0:11 



0:12 



(61) 



which is consistent with eas. (|56p . (|59[) . In particular, we 
have: 

4 f°° 

K,n + K,22 = ~Y k2 J d 1 p M- (62) 
Moreover, at high k we have the asymptotic behaviors: 

dqP L0 (q) (63) 

2 4tt 



k — > 00 : S 0;11 ~ £q;21 



34n r°° 
5 3 . 



^0;22 



^0:12 



5 3 



dqP La (q) (64) 



Of course, we can check that this is consistent with ea. (|62l) . 

We display in Fig. Q] the self-energy term £g (fc) as a 
function of wavenumber k, from which E can be obtained 
using eqs. ([59")) . (fBTj) . All components Ej^ ■ are of the same 
magnitude and negative except for Sq-ii which is positive 
at large scales k 0.2h Mpc -1 . All terms grow as k 2 
at high wavenumbers and we recover the asymptotic be- 
haviors (|63|) - (f64|) . Note that the self-energy decreases at 
large scales so that the two-point functions obtained from 
cqs.(j42 |) - (|44l) will converge to the linear asymptotics 

We show in Fig. [2] the evolution forward over time rji 
of the self-energy £0 (A;; 771, 772) for time 772 = —1.4 (i.e. 
z 2 = 3) and wavenumbers k = 0.1, 1 and 10 x h Mpc -1 
(from bottom to top). The absolute value of the self-energy 
is larger for higher k in agreement with Fig. [TJ All com- 
ponents are negative (except for S 0: n at r/i > —0.8 in the 
case k = O.lh Mpc -1 ) and exhibit a smooth dependence 
with time following eq. (|5T)|) which goes to e 2ri1 at large 771 . 

We show in Fig. [3] the evolution backward over time 
772 of the self-energy Eq (A;; 771, 772) for time 771 = (i.e. 
z\ = 0) and wavenumbers k = 0.1, 1 and 10 x h Mpc -1 
(from bottom to top). The behavior agrees with Figs. ll|2l 
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Fig. 2. The self-energy £o(fc; 771, 772) as a function of for- 
ward time 771, for 772 = —1.4 (i.e. z-i = 3) and wavenum- 
bers k = 0.1, 1 and 10 x h Mpc -1 . The line styles are as 
in Fig. [TJ The diagonal components vanish for 771 < 772 
whereas the off-diagonal components vanish for 771 < 772. 
All terms are negative except for £o;ii which becomes 
positive at 771 > —0.8 for k = O.lh Mpc -1 . 
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Fig. 3. The self-energy T, (k; 771, 772) as a function of back- 
ward time 772, for 771 = (i.e. z\ — 0) and wavenumbers 
k = 0.1, 1 and 10 x h Mpc -1 . The line styles are as in 
Fig. Q] All terms are negative except for £o ; ii which be- 
comes positive at 772 < —0.6 for k = O.lh Mpc -1 . 



Following eq. (|59"l) So converges to a constant for 772 — > 
—00. 

Finally, using the factorized form (|21| of the auxiliary 
two-point correlation the self-energy Hq can be written as: 



IIo (si, a: 2 ) = Mki + k 2 )e 2 " 1+2 ^no(fci), 



Fig. 4. The self-energy Ho(k) as a function of wavenumber 
k. At high k all components are very close and positive. 
At k < 0.06ft. Mpc -1 the off-diagonal components become 
negative. 



with: 

n : 4l i 2 (fc) 



jihhh 



,(-q.-k 



J d q7u yiJ2 (q,k-q) 

Fq)P M (g)P M (|k-q|). 



(66) 



x 7 l2 ;i 1 ,; 2 * 

In agreement with the symmetry f|40[) we can check that 

no;2i(A) = n 0; i 2 (fc). 

Since the time-dependence (|65j) is quite simple we only 
display the self-energy term U (k) of eq. l[66|) as a function 
of wavenumber k in Fig. |U At high k all components are 
very close and positive. The dependence on wavenumber k 
follows the smooth growth at small scales of the logarith- 
mic power A| /0 (fc) = fc 3 Pio(fc). As for £0 the self-energy 
term IIo decreases very fast at low k so that the two-point 
functions obtained from eqs. (|42|) - (|44|) will converge to the 
linear asymptotics at large scales. 

5.2. Response function 

From the self-energy obtained in sect. I5.11 we compute the 
response function from eq. (|43|) . Thanks to the Heaviside 
factor within the term Ho(x, z) the r.h.s. only involves 
earlier times hence eq. (|43p can be integrated forward over 
time 771 at fixed 772 to give R(k; 771, 772). Note that each 
wavenumber k evolves independently thanks to the Dirac 
factors <$D(ki — k 2 ) within both £0 and R as in eq. ([50|) . 
Thus all processes associated with mode-coupling between 
different wavenumbers k are contained in the calculation 
of the self-energy £ (and n ), as in eqs. (|60"|) . ((66)) . 
Eq. ([43| also writes for 771 > 772: 

m „ „ n , „ (67) 



dR 2 



( 65 ) drn 



dm 



— i?2 — £():11--Rl + £o;12--R2 



1 



-i?2 — £():21--Rl + £o;22--R2 



(68) 
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where (i?i,i? 2 ) = (i?n,i?2i) or (i?i 2 ,i? 22 ), indeed in this 
case the two columns of the matrix R are not coupled. 
Using eqs. fST))) and (foTj) the first eq. (|57|) also reads: 



dRi 

dm 



(571,172) - #2(771,772) = 



d V |(e 2t »£+ u +e-" 1 / 2+5 "/ 2 S+ 22 )i? 1 ( ?7 , 772 ) 
f (e 2,,1 S+ 12 - e-^+^X+^R^m)} (69) 



Taking advantage of the simple dependence on time 771 of 
the r.h.s. it is possible to eliminate the integrals over 77 by 
taking two derivatives of eq. (|69[) with respect to 7/1 . This 
yields: 



d 3 R! 
drj\ 

2m 



Z d 2 R 1 
2 9?7 2 

(^o-ii 



<9 2 i? 2 ^ 3 <9i? 2 

drjl drj! 



^22) 



dR ± 



2 am 



i?2 = 



a?7i 2 

In a similar fashion, one obtains from eq. (|68p : 
d 3 R 2 3 9 2 i?i d 2 R 2 9 9i?i 7 dR 



E 012 R 2 



(70) 



9»7i 



,2»7i 



2 <9?/ 2 

(^0:11 



9?7 2 4 dr]i 
,+ ,dR 2 



4 % 



— Ri — — i? 2 — 
2 2 



J 0:22 



(71) 



These equations apply to r/i > 7y 2 and need to be supple- 
mented by initial conditions obtained from eas. (fB7 |) -(f55 |) 
and eq. (|23[) which yield at point 771 = 77+ : 



X 1 



r> r> d-Bi dR 2 



1 § 



=2*71 



(^0:11 



^0-22) 



(72) 
(73) 



for the pair (i?n,i?2i 
l r = (0 11-i- 



and: 



+ e 2 ^(E+ u + S+ 22 )) (74) 



for the pair (i?i 2 , i? 22 ). One can easily check that the lin- 
ear response Rl of eq. (|2"4"|) is a solution of eqs. (|7D|) - ([7T|) 
with Eq = 0. The advantage of these differential equations 
over the integro-differential eqs. ([6T)) - (p5|) is that they no 
longer require the computation of integrals over past val- 
ues of the response R. Therefore, one does not need to 
save the values of R on a grid with a very fine mesh and 
one can advance over large time steps 771. Moreover, each 
time-step runs faster since one does not need to compute 
lengthy integrals as in eqs. ([6T|) - ([rJ5|) . We use an adapta- 
tive Runge-Kutta algorithm to solve numerically between 
two grid points the eqs. ([70)l - (|7T|) , written as a system of 
coupled first-order linear equations over the components 
of the vector X of eq. (|72"]) . 

We can note that eqs. ([7D)) -(|7i p are somewhat similar 
to Bessel functions in terms of the scale factor a\ = e'' 1 
but are of degree three rather than two. In the small-scale 
limit I £ q I — > 00 we can look for an asymptotic solution of 
the form: 



Ri.Avi) = PiaimV' 



(75) 



z 2 = 3, 77 



1.4 




Fig. 5. The non-linear response R(k; rji, ?y 2 ) as a func- 
tion of forward time 771, for 772 = —1.4 (i.e. z 2 = 3) and 
wavenumbers k = 1 (left panel) and 10 x h Mpc -1 (right 
panel). We also plot the linear response Rl which shows 
a simple exponential growth. 



where the functions Pi(rji) may be expanded as: 

= ^ + + + (76) 
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Fig. 6. The non-linear response R(k\ 771, ?7 2 ) as a function 
of backward time ?y 2 , for 771 = (i.e. z\ = 0) and wavenum- 
bers k = 1 (left panel) and 10 x h Mpc -1 (right panel). 
We also plot the linear response Rl which shows a simple 
exponential growth. 
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and i 
into eas.([70|)-([7l 

UJ = ± 



|Sq |. We obtain at order ui 3 after substituting 



4tt 



EJu - K,22 = ±k ( y / dgP i0 (g) ) (77) 



1/2 



where we used eq.(|62j). At order u we obtain: 



dpi 



(0) 



(o) 

P 2 



9p 



(0) 
2 



3 (o) 



(0) 



= 







(78) 



(79) 



where we used the asymptotic behavior ([53")) - ([M]) to sim- 
plify factors of the form Sj^/w 2 . By comparison with 
eas. (|67l) - (|68)) . we see that these equations are precisely 
the equations obeyed by the linear response Rl = Ro- 
Moreover, the initial conditions at 771 = 772 are precisely 
the same if we redefine at order 0: 

R^(rn) = pf ) (7 7l )cos[ W (e" 1 -e 1 ' 2 )] 

= Pi°\ax) cos[o;(ai - a 2 )]. (80) 

Therefore, at high k the "envelope" pf^i^i) of the non- 
linear response i?y (771) follows exactly the linear response 
Rl. In addition, the response R exhibits periodic oscil- 
lations over the time variable a\ with a frequency which 
increases at high wavenumbers as to oc k. Therefore, at 
one-loop order within the direct steepest-descent expan- 
sion the memory of initial conditions is not really erased 
at small scales or late times since the amplitude of the 
response R does not decay but only exhibits strong oscil- 
lations. 
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Fig. 7. The non-linear response function R(k; 771, 772) as a 
function of wavenumber k, at times 771 = 0, 772 = —1.4, i.e. 
Z\ — 0, Z2 = 3. The horizontal lines are the linear response 
Rl of eq. which is independent of k 

We first display in Fig. [5] the evolution forward 
over time 771 of the non-linear response R(k; 771, 772) for 



wavenumbers k — 1 (left panel) and 10 x h Mpc -1 (right 
panel). We do not show the case k = O.lh Mpc -1 where 
the response R closely follows the linear response Rl 
given in eq . (|24|) . In agreement with the analysis performed 
through eqs. ([75l) - ([80|) we find that the non-linear response 
R exhibits oscillations ~ cos(wai) with a frequency lo 
which grows at higher k and an amplitude which follows 
the linear response Rl- Moreover, all components Rij os- 
cillate in phase. 

Next, we show in Fig. [S] the evolution backward over 
time 772 of the non-linear response R. It shows a behav- 
ior similar to the forward evolution displayed in Fig. [5] 
with respect to 771. At large scales we recover the smooth 
increase (in absolute values) at early times 772 of linear 
theory (|24|) whereas at small scales we obtain fast oscilla- 
tions ~ cos(wa2) related to eqs. (|T5)) - (15Ul) with an ampli- 
tude which follows again the linear response Rl- 

Finally, since at equal times the response function 
obeys eq. (|23|) we only display in Fig. [7] the response at 
unequal times (771 = 0,772 = —1.4). At low k we recover 
the linear response Rl- At higher k above 0.2h Mpc^ 1 the 
non-linear response obtained at one-loop order within this 
steepest-descent approach departs from the linear predic- 
tion and shows oscillations with increasingly high frequen- 
cies while their amplitude follows the linear response Rl 
at high k. This behavior is a result of the oscillatory be- 
havior with time 771 displayed in Fig. [5] above and analyzed 
in eqs.([75|-(j80l). 



5.3. Two-point correlation 

Finally, from eq. (|51j) we compute the non-linear two-point 
correlation function G obtained through the steepest- 
descent method at one-loop order. We compare our re- 
sults in Figs. rSlfTTIat redshifts z = 0, 3 with a fit G n b from 
numerical simulations (Smith et al. 2003), the linear pre- 
diction Gl feq. ((2Tj) ') and the usual one-loop result Gii oop 
obtained from the standard perturbative expansion over 
powers of the linear density field. The standard one-loop 
power-spectrum can be written (e.g., Jain & Bcrtschinger 
1994, Scoccimarro & Frieman 1996a, b) as: 



Pi 



Hoop 



= Pl 



P: 



22 



A3, 



(81) 



where Pl is the linear contribution and P22 and P13 are 



of order P| with: 



p22(fc; 771,772) 
and: 



e 2„ 1+2 „ 2 y dqPio((? )p LO ( k _ q ) 
x2F 2 (s) (q,k-q) 2 , (82) 



-(e 3 " 1+ " 2 + e " 1+3??2 )P L o(fc) 
x J dqP i0 (<7)6F 3 W (q,-q,k). 



(83) 
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The angular kernels F 2 
couplings a and (3 of eq 
et al. 1986): 



is) 



(s) 

and F£ ' are obtained from the 
i}. In particular we have (Goroff 



^ 2 (s) (ki,k 2 ) 



ki-k 2 

2kik 2 



k2 
fcl 



ki-k 2 

k\k 2 



.(84) 



On the other hand, after integration over angles in eq. 
one obtains (Makino et al. 1992): 



Pi3(k]r]i,r) 2 ) 



1 



= 3l71+7/2 



= 171 +3^2^ 



Pw(k) 



xvr J dqP L0 {q)q 2 

(g 2 -fc 2 ) 3 (7g 2 + 2fc 2 ) 
+ 42k 3 q 5 



6fc 6 - 79fcV + 50&V - 21q 6 



In 



63k 2 q 4 
k + q 



(85) 



From eqs. (1821) - ([85)1 we can see that at large wavenumbers 
k — > 00 the integrals over q in eqs. (|52")) . (|85[) are domi- 
nated by the region where the slope n of the linear power- 
spectrum is of order n ~ — 1. This yields: 

fc -» 00 : P lloop ~ -e 2 ^+2m ( cos h( m - m ) - 1) 

xk 2 P L0 (k)^- J dqP La (q). (86) 

Therefore, at unequal times 771 ^ 772 the standard one-loop 
prediction becomes negative at high k and the logarithmic 
power (as defined below in eq. ([87| ) grows as A 2 loop oc 
— e 3m k 2 with k and 771, in agreement with Figs. [SJ [TT] 
Note that the leading term (|86[) vanishes at equal times 
771 = r]2 ■ This is related to the cancellation of the infrared 
divergences which appear for linear power-spectra Pj, oc 
k n with n < — 1 when one considers equal-times statistics 
(Vishniac 1983; Jain & Bertschinger 1996). 

The dependence e 3m k 2 merely reflects the order P 2 
of the perturbative expansion. In order to go deeper into 
the non-linear regime one may consider higher order terms 
but the latter grow increasingly fast and one would need at 
least partial resummations to obtain a well-behaved series 
(e.g. Crocce & Scoccimarro 2006a). This is precisely what 
the large- TV expansions discussed in this article attempt 
to perform. 

We first show in Fig. [8] the evolution over time 771 of 
the density two-point correlation Gn(fc; 771, 772), which is 
symmetric with respect to 771 <-> 772 . We actually display 
the "logarithmic power" A 2 (fc; 771, 772) defined as in eq. (fT9|) 
by: 



A 2 (fc; 771,772) =4ttA: 3 Gii(A:;77i, 772). 



(87) 



At equal times 771 = 772 it is equal to the usual logarith- 
mic power and we have A 2 > 0. However, at different 
times 771 7^ 772 we can have Gn < whence A 2 < 0. 
We present the cases k — 1 (left panel) and k = 10 x h 
Mpc -1 (right panel) in Fig. [5] for 772 = —1.4. For clarity 
we only display the density-density correlation G\\. Other 
components G^ show a similar behavior. 

The oscillations show that Gn can indeed be nega- 
tive at unequal times. At these small scales the non-linear 
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Fig. 8. The density two-point correlation Gn(fc; 771, 772) as 
a function of time r)\, for 772 = —1.4 (i.e. z 2 = 3) and 
wavenumbers k = 1 (left panel) and 10 x h Mpc -1 (right 
panel). We plot the logarithmic power A 2 = 47rfc 3 Gn. We 
display the linear power A| of eq. (|2"Tj) (dashed line), the 
usual one- loop perturbative result A 2 loop of eq. (|8Tj) (dot- 
ted line), the full non- linear power A 2 from eq. ([5T| (solid 
line) and the contribution A 2 from eq . (f52|) (dot-dashed 
line). Note that Gn whence A 2 can be negative, as shown 
by the oscillations displayed in the figure. 



power A 2 is already significantly different from the lin- 
ear power A| but we can check in the left panel that 
A 2 matches the usual one- loop power A 2 loop of eq. (f8Tj) at 
small A 2 . As recalled in sect. !4.3l the one-loop large- N ex- 
pansions indeed match the usual perturbative expansion 
up to one-loop order (i.e. Pi). At later times while A 2 loop 
grows very fast as e 3m from cq. ([83|) the prediction A 2 
of the one-loop steepest-descent approach remains well- 
behaved of order A 2 L . Thus, we can see in the right panel 
that at very small scales where A 2 loop has become exceed- 
ingly large A 2 remains well-controlled. In fact, following 
the behavior of the response function analyzed in sect. 15.21 
it exhibits fast oscillations with an amplitude close to A 2 . 
We also show in Fig. [H] the power A 2 obtained from the 
first term alone of eq. (p)Tj) which is detailed in eq. (|52"|) . This 
corresponds to the "non-linear transport" of the initial lin- 
ear fluctuations without taking into account the fluctua- 
tions II generated at all times by the non- linear dynamics. 
We can see that this term is sub-dominant in the mildly 
non-linear regime (left panel) but happens to be again of 
same order as A 2 in the highly non-linear regime. 

Next, we show in Fig. [9] the logarithmic power A 2 (fc) 
today at equal redshifts z± — z 2 — (i.e. equal times 
771 =772 = 0). In agreement with eqs.([5"l" j) -(f52" |) and the 
subsequent discussion both A 2 and A 2 are positive (A 2 is 
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Fig. 9. The logarithmic power A 2 (A;) of eq. (fT9|) at redshift 
z = 0, that is for equal times 771 = 772 — (zi = Z2 = 0). 
We display the linear power A 2 (dashed line), the usual 
one-loop perturbative result A 2 loop (dotted line), the full 
non- linear power A 2 from eq. (p)Tj) (solid line) and the con- 
tribution A 2 from eq. ([52")) (dot-dashed line). We also show 
the fit A 2 b ( upper solid line) from numerical simulations 
(Smith et al. 2003). All quantities (including A 2 which is 
the square of an oscillating function) are positive. 
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Fig. 10. The logarithmic power A 2 (fc) as in Fig.[9]but at 

redshift z = 3. 

the square of an oscillating function). At large scales A 2 
converges to the linear regime A^. At smaller scales we 
can check again that A 2 first follows the usual one-loop 
power A 2 loop and next goes back to A| whereas A 2 loop 
keeps increasing very fast in magnitude. Unfortunately, 
A 2 does not agree with the fit A 2 b from numerical simu- 
lations (Smith et al. 2003) on a larger range of k as com- 
pared with A 2 loo (although the rough agreement between 



Fig. 11. The logarithmic power A 2 (fc) as in Fig. [9]but at 
unequal times (771 = 0,772 = —1-4) (i.e. z\ = 0, z 2 = 3). 

A 2 b and A 2 loop at small scale is purely coincidental). On 
the other hand, whereas adding higher-order terms in the 
standard perturbative expansion may not improve much 
the agreement with N-body simulations (especially since 
it would lead to increasingly steep terms at high k) the 
series obtained in the large-^V approach is likely to be 
well-behaved as various terms do not "explode" at small 
scales. However, this would require intricate calculations. 
We can note that although A 2 is of order A 2 one can- 
not neglect the fluctuations II generated by the non-linear 
dynamics and the sum yields a smooth power A 2 . 

We also display in Fig.Qjllthe logarithmic power A 2 (fc) 
at equal redshifts z\ = Z2 = 3. We recover the same be- 
haviors as those obtained at redshift z = 0. 

Finally, we show in Fig. QT] the logarithmic power 
A 2 (fc; 771, 772) at unequal times 771 =0,772 = —1-4, that is 
at redshifts Z\ = 0, = 3. We find again that A 2 matches 
the usual one-loop power A 2 loop at large scales and follows 
its change of sign at k ~ O.8/1 Mpc -1 . Note indeed that at 
unequal times A 2 need not be positive. At smaller scales, 
the steepest-descent large-A result A 2 departs from the 
fast growing A 2 loop and shows a series of oscillations of 
moderate amplitude, which again follow A 2 . These fea- 
tures suggest that the first change of sign in Fig. [TTJ which 
is common to both A 2 ^ and A 2 loop , may be real and not a 
mere artefact of the one-loop expansion. Note that in spite 
of the oscillations with time seen in Figs. I8|lll the large- 
N approach automatically ensures that A 2 > at equal 
times thanks to the structure of eq. (|5ip . This property 
holds at all orders over 1/N since ea.([ST|) is independent 
of the expressions of the self-energies S, II. 

6. 2PI effective action approach 

We now present the results obtained at one-loop order 
from the 2PI effective action approach recalled in sect. 14.21 
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Thus, we need to solve the system of coupled equations 
(|42] ) -(|36 |) where in the expression Ipp ft -p)] ) of the self- 
energy the auxiliary two-point functions Go and R are 
replaced by G and R. Thanks to causality, which leads to 
the Heaviside factor 9(rii — rj 2 ) of eq. (|5D|) within both R 
and S, we solve the system (l4"2j) - (l4T>l) by moving forward 
over time. 

Thus, we use a grid rf n > with n = 1, .., N v for the time 
variables. At the earliest time-step r/ 1 ) we initialize all ma- 
trices by their linear value at 771 = 772 = V ■ Next, once we 
have obtained all two-point functions up to time rj^ (i.e. 
over 771 < rf n > and 772 < T)( n ', initially n = 1) we advance 
to the next time-step 77^™ +1 ^ as follows. The response R at 
equal times 77! = 772 = 77(" +1 ) is first obtained from eq. p5|) . 
Next, we move backward over time 772 = r]( n \ .., T)^ 1 ' at 
fixed 771 = 77(" +1 ) by using eq. (|l3")) and eq. (|^5|) . That is, 
to compute R(r]^ n+1 \r]^) with i = n, .., 1 we use for each 
value of i the integro-differential eq.(!43|) to move over 771 
from 771 = 77'™) up to 771 = 77( n+1 ) at fixed 772 = 77W. 
Thanks to the Heaviside factors the r.h.s. of eq. (|43|) only 
involves 2(771,77) and ^(77,772) with 771 > 77 > 772 which 
are already known. Besides, once 

j R( r? («+i) ;?7 W) has been 
obtained we compute X(t/™ +1 ) , 77W) from eq. (|45f before 
moving downward to step i — 1. Since R and S vanish for 
772 > 771 we have actually obtained in this fashion R and 
£ for all 771,772 < ?7 ( " +1) (indeed R{r) (i) , ^ n+1 ^) = for 
i < n). 

Next, we compute II and G from eq. (|46|) and eq. ([5Tj) . 
At fixed 771 = rj^ n+1 ^ we now move forward over time 



where <5y is the Kronecker symbol. The basis functions 
X^ n ' (k) define our interpolation scheme. We could use ba- 
sis functions which all extend over the full range < 
k < 00 (i.e. a spectral method built from cardinal func- 
tions) but we prefer here to use simple finite elements 
with x {n) {k) = if k < fc^" 1 ) or k > /c( n+1 ). For conve- 
nience we actually use triangular functions over the vari- 
able k = (k - l)/(k + 1) which maps ]0, oo[ to ] — 1, 1[: 



Jn-l) 



X W - S K (n+1) _ K 



K (n+l) _ K (n) 



if K^-V < K < K (n) 
if «W < K < 



(90) 



which is equivalent to use linear interpolation over k for 
G(k). Then, we can take advantage of the fact that the 
kernels r y-. ii i2 (ki,k 2 ) in eqs. ([53|) - ([54|) do not depend on 
time to compute once for all the integrals over wavenum- 
ber needed for the self-energy £ and II. Thus, we define 
the quantities: 

S^'" 2) = 4 J dq 7 ^(q,kW-q)7^(kM q-kM) 

xx (ni) (9)x (na) (|k (n) -q|) (91) 

and: 

n^ 1 '" 2 ' = 2 J dq 7 « 1 (q,kW-q)7^(-q ) -kW+q) 

xx (ni) («)x ( " 2) (|k (n) -q|) (92) 
where we noted v = (i; 11,72) the indices of the vertices 



m 



,77 (n) . Indeed, the r.h.s of eq.jH]) only in- Ji;i lt i 2 - These quantities £ 



(n;«i, n 2 ) 



and II 



(n;ni,n 2 ) 



only de- 



volves 11(77^, rj' 2 ) with 77J < 771 and r)' 2 < 772 which is already 
known. Besides, once G(77(™ +1 ) , 77W) has been obtained we 
compute n(?/ rl+1 ), 77W) from eq. (|46l) before moving for- 
ward to step i + 1 and we use the symmetry (|49|) to derive 
G^W^^+i)) as well as n^M,??^ 1 )). 

In fact, at each step i for (?/" +1 ), 77W) the integrals 
in the r.h.s. of eqs. ([43|) . (|51[) also involve the values 
R(r]( n+1 \ri^) and G(rj( n+1 \ 77W) which are being com- 
puted if we use the boundary points in the numerical 
computation of the time- integrals. Rather than using open 
formulae for the integrals we perform a few iterations at 
each time-step rf n ^ . This procedure converges over a few 
loops. 

Finally, to speed-up the numerical computation we use 
finite elements over Fourier space k to store the structure 
of the self-energy terms £ and IT Thus, since all matrices 
only depend on the wavenumber modulus k = |k| we use 
a grid few with n = 1, .., Nk- Next, in order to interpolate 
for all values of k we can write for instance any matrix 
such as G{k) as: 



pend on the vertices 7 s given in eqs. (fll)|) - (ri"tj|) and on the 
grid which we use over wavenumbers hence they are com- 
puted at the beginning of the code within initialization 
subroutines. Next, the self-energy at a given time 771,772 is 
obtained as: 

E ili2 (fcW 5771 ,772) 



E 



(93) 



with vi = (ii;j'i,j2) and v 2 — (h;i2,h), and: 



xGJ 2 2 ) (77 1 , 772) 



n=l 



with: 



G (n) = G ( fc («)) and = §, 



131 



(89) 



(94) 

with v\ = (ix;ji,jz) and v 2 — faihth)- Therefore, us- 
ing eqs.([93 ]) - ([94l) we obtain the self-energy by advancing 
over times 771 , 772 as described above without performing 
new integrations over wavenumbers at each step (they are 
replaced by the discrete sums over 77,1,712 with weights 

v (n:ni,n 2 ) TT (n:"l,n2)\ 



6.1. Response R and self-energy S 

We first display in Fig.[12]the evolution forward over time 
771 of the response R(k; r\\, 772). We can see that the non- 
linear response exhibits oscillations as for the steepest- 
descent result of Fig. [5] but its amplitude now decays at 
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z 2 = 3, 77 2 =- 1.4 




Fig. 12. The non-linear response 771 , 772) as a func- 
tion of forward time 771 , for 772 = —1.4 (i.e. z 2 — 3) and 
wavenumbers k = 1 (left panel) and 10 x /i Mpc -1 (right 
panel), as in Fig.[5]but for the 2PI effective action method. 
For comparison we also plot the first half-oscillation of 
the response obtained in sect. [5] from the direct steepest- 
descent method (curves labeled "sd" ) which was shown in 
Fig. El 



large times 771 instead of following the linear envelope. This 
can be understood analytically from the following simple 
model. Let us consider the equation: 



OR 

dm 



dr? e^+^i, 77)^(77, 772) 



(95) 



'/2 



for a one-component response R^x,^). This is a simpli- 
fied form of the system (|4"3"]) . (jl5)) where we used the linear 
time-dependence (f2~Tj) for G to obtain a closed equation for 
R. The parameter a < represents the amplitude of the 
self-energy £ at the wavenumber of interest. Then, the 
linear regime Rl corresponds to a = and yields the 
solution: 



m>V2- Rl(vi,V2) = ^, 



(96) 



where we used the initial condition R(j], 77) = 1 at equal 
times. Next, the steepest-descent method presented in 
sect. 14.11 and sect. [5] amounts to replace R(r]i 7 r]) by 
Rl{i1i>v) m to the integral in the r.h.s. of cq. ([95)) . This 
yields a linear equation for R which can be transformed 
into a differential equation as in sect. 15.21 by taking the 
derivative of eq. (|95p with respect to 771 . Since we only have 
one component one differentiation is sufficient to remove 
the integral and we obtain: 



d 2 R 
dij 2 



^ = ae^R 



whence 



d 2 R 



= aR, 



(97) 



where we introduced the scale-factor a\ — e rn as the new 
time- variable. Therefore, we obtain the solution: 



R(ai, 02) = cos[w(ai — 02)] with uj = \/—a. 



(98) 



Thus, as in sect. 15.21 we obtain within the steepest-descent 
method an oscillating response function with an amplitude 
given by the linear response Rl and a frequency over the 
time-variable a\ given by \J —a. Finally, in the 2PI effec- 
tive action method we need to solve the non-linear equa- 
tion which also reads: 



OR 



a / do R{a\, a)R(a, a 2 ) 



(99) 



where we changed variables from 77 to a — e ri . We can look 
for a solution of eq. ([99|) of the form R(a\, 02) = R(ai — 02) 
so that the r.h.s. of eq. ([99|) becomes a simple convolution: 



dR 
da 



da' R(a - a')R(a'). 



(100) 



Introducing the Laplace transform R(s) defined by: 

R(s) = da e~ sa R(a), (101) 
Jo 

and taking the Laplace transform of ea. (|100p we obtain: 

sR{s) - i?(0) = aR(s) 2 . (102) 

Using R(0) = 1 and requiring that R(s) vanishes for s — > 
+00 yields: 



R(s) 



2uj 2 



(103) 



where we introduced u> as defined in eq. (|98p . This is a 
well-known Laplace transform and eq. (|103p directly gives: 



R(a 1 ,a 2 ) 



Ji[2oj(ai — 02)] 
uo(ai ~ a 2 ) 



(104) 



where Ji is a Bessel function of the first kind. Thus, we see 
that within the 2PI effective action method the response 
function exhibits oscillations as in the direct steepest- 
descent method. However, their amplitude is no longer 
given by the linear response and becomes much smaller at 
late times along with a higher frequency. For the simple 
model (|95j) the amplitude of the response function decays 



^ /2 

as a 1 instead of being constant while the frequency be- 
comes twice larger. This explains the qualitative behavior 
seen in Fig. 1121 In particular, note that in agreement with 
the simple model (l95l) - (|104p the comparison of the right 
panel of Fig. [T3] with Fig. [5] shows that the frequency of 
the oscillations is indeed about twice larger in the non- 
linear regime for the 2PI effective action result than for the 
steepest-descent prediction. Of course, the result shown in 
Fig. [12] is more complicated than the simple model (|95[) 
since the two-point correlation G which enters £ whence 
R also depends non- linearly on R through eq. (!42p or the 
more explicit expression (|5ip . However, as we shall check 
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Fig. 13. The non-linear response i?(/c; rji, 772) as a function 
of backward time 772, for 771 — (i.e. Z\ — 0) and wavenum- 
bers k = 1 (left panel) and 10 x h Mpc" 1 (right panel). 
We also plot in the left panel the first half-oscillation of 
the response obtained in sect. [5] from the direct steepest- 
descent method (curves labeled "sd" ) which was shown in 
Fig. El 
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Fig. 14. The non-linear response function R(k; 771, 772) as a 
function of wavenumber k, at times 771 = 0, 772 = —1-4. We 
also plot the first half-oscillation of the response obtained 
in sect. [5] from the direct steepest-descent method (curves 
labeled "sd" ) which was shown in Fig. [7] 



in sect. 17.31 the dependence of S on G does not change 
much the behavior of the response R (see Fig. |2"5]) . 

On the other hand, we can see in Fig. [12] that the first 
quarter of oscillation, from 1 to 0, of the response R fol- 



lows rather closely the curves labeled "sd" obtained from 
the steepest-descent method presented in sect. [5] which 
were also shown in Fig.O Besides this agreement remains 
valid at quite large k (right panel) where two-point func- 
tions obtained from both methods are generically very dif- 
ferent as shown by the figures below. This property can 
be understood from eq. J56J) which shows that the depen- 
dence on wavenumber £(fc;77, 77) oc k 2 of the self-energy 
£ at equal times is identical for both approaches (at one- 
loop order). Moreover, eq. ([56|) shows that the normaliza- 
tion is governed by the value of the two-point correlation 
G22(fc'; Vi v) a t the point where the logarithmic slope is 
n — — 1 (this integral actually corresponds to the mean 
square velocity (v 2 )). At redshift z = this wavenum- 
ber is still within the linear regime (for the SCDM cos- 
mology which we consider here) therefore the self-energy 
£ at equal times r\\ = 772 obtained within the steepest- 
descent method and the 2PI effective action approach are 
very close until z = 0. Then, the early time-evolutions at 
771 ~ 772 of the response R obtained within both methods 
from eq. (|43"]) are very close. We can see from Fig. [12] that 
this agreement holds until the response R first vanishes. 
Beyond this point the steepest-descent method yields in- 
creasingly large oscillations (Fig. [5]) whereas the 2PI ef- 
fective action yields small oscillations which decay to zero 

(Fig. ED. 

Next, we show in Fig. [T3] the evolution backward over 
time 772 of the response R. In a manner consistent with 
Fig. [12] it exhibits oscillations which are strongly damped 
at large time separations |ai — 02!, whence at early times 
772 , while the frequency is also larger than for the steepest- 
descent result shown in Fig. [6] Again the behavior at 
nearly equal times 772 — 771 is close to the prediction of 
the steepest-descent method presented in sect. [5] until the 
response first vanishes. The response at unequal times 
(771 = 0, 772 = —1.4) is shown as a function of wavenumber 
k in Fig. [14] At large scales we recover the linear regime 
whereas at small scales we obtain oscillations which show 
a fast decay into the non-linear regime. This is consistent 
with the time-dependence displayed in Figs. [T2] [13] 

Therefore, in contrast with the steepest-descent ap- 
proach we find that within the 2PI effective action method 
the memory of initial conditions is in some sense erased 
as the response function decays for large time separation. 
This agrees with the intuitive expectation that within the 
real collisionless gravitational dynamics the details of the 
initial conditions are erased at small scales. Indeed, after 
shell-crossing one can expect for instance that virializa- 
tion processes build halos which mainly depend on a few 
integrated quantities which characterize the collapsed re- 
gion (such as the mass, initial overdensity, angular mo- 
mentum,..) and exhibit relaxed cores (e.g. isotropic veloc- 
ity distributions), as suggested by numerical simulations. 
However, we must note that the system studied in this pa- 
per is defined from the hydrodynamical equations of mo- 
tion (H])-© which break down at shell-crossing. Therefore, 
there is no guarantee a priori that such small-scale relax- 
ation processes should come out from eqs.([T])-([3]). It ap- 
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pears through the 2PI effective action method presented in 
this paper and Figs. [T2]fT4l that this is actually the case. In 
fact, studying the large- fc behavior of the kernel K s Crocce 
& Scoccimarro (2006b) managed to resum all dominant di- 
agrams in this large-fc limit and obtained a Gaussian decay 
R(k) ~ e~ k . In our case, the simple model ([95]) suggests 
that the one-loop 2PI effective action method only yields 
a power-law decay as in ea. (|104l) . Since at the one-loop 
order considered in this article we do not resum all dia- 
grams it is not really surprising that we do not recover the 
"exact" Gaussian decay. However, note that the large-TV 
expansion does not give a mere Taylor series expansion of 
such a Gaussian decay (which would blow up at large k 
or large times) but already captures at a qualitative level 
the damping of the response R at small scale or at large 
time separation thanks to the non-linearity of the evolu- 
tion equation obeyed by R. 



z 2 =3, 77 2 = - 1.4 
10 1 + Zj l 10 1 + Zj l 



er 10 











// 






k=10 h 


/Mp$/' 

/// / 


: k=l h 


/Mpc 




f 


/sd : 




sd/; 

//, 




















/ ..• \\\- 


I 1 I 1 1 1 







-2 -1 0-2-10 

Fig. 15. The self-energy E(fe; 771, 772) as a function of for- 
ward time 771, for 772 = —1.4 (i.e. Z2 = 3) and wavenum- 
bers k — 1 (left panel) and 10 x h Mpc -1 (right panel). 
The curves labeled "sd" are the prediction of the steepest- 
descent method as in Fig. [2] 

We display in Fig. [15] and Fig. [16] the evolution with 
times r]i and 772 of the self-energy E. At large scales 
(k < O.lh Mpc -1 ) and small time separations we recover 
the results obtained within the direct steepest-descent 
method but at smaller scales we obtain a series of oscilla- 
tions with an amplitude which decreases at large time sep- 
arations. This follows from the behavior of the non- linear 
response R analyzed in Figs. [T2iri4l above. Indeed, the self- 
energy E depends linearly on the response R, see eq. (f93]) . 
Again the behavior at nearly equal times 771 ~ 772 is sim- 
ilar for both approaches (see in particular the left panel 
of Fig. [15]). The right panel of Fig. [T5l shows that at high 
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Fig. 16. The self-energy E(fc; 771, 772) as a function of back- 
ward time 772, for 771 = (i.e. z\ = 0) and wavenum- 
bers k = 1 (left panel) and 10 x h Mpc -1 (right panel). 
The curves labeled "sd" are the prediction of the steepest- 
descent method as in Fig. [3] 
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Fig. 17. The self-energy E(fe) as a function of wavenumber 
k at unequal redshifts z\ = 0, zi = 3. The curves labeled 
"sd" are the prediction of the steepest-descent method. 

wavenumbers the envelope of the oscillations of E decays 
as a power-law of a\ = e Vl at late times. This agrees with 
the analysis of the simple model ([^5]) . In a similar fashion, 
Fig. ll7l which displays the dependence on wavenumber k of 
E(fc; 771, 772), at fixed times 771 = 0, 772 = —1-4, matches the 
results obtained from the direct steepest-descent method 
at low k (except for E12) and exhibits decaying oscilla- 
tions at high k. We can note in Fig. [17] that the term E12 
(dotted line) does not match the steepest-descent result 
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in the limit k — > 0. This is due to the behavior of the asso- 
ciated kernels E^™;,™ 1 '™ 2 ^ of eq. ([9Tj) . More precisely, using 
eqs. ([T5j) - (fr5|) . eq. (f53"l) reads in this case: 

[k.(k-q)](k.q) 



z 2 = 3, 7} 2 — — 1 A 



10 1 + z j 1 10 l+z 1 i 



10 



Eia(fc) = / dq i? n (<z)G 21 (|k-q|) 



-#i 2 (g)G 22 (|k-q|) 



/c 2 |k-q| 2 
k.(k-q)] 2 g 2 
/c 2 |k-q| 4 



+i22i(g)Gu(|k-q|) 



(k-q) 2 

k 2 q 2 



^ 1 

CM 
< 



0.1 - 



Er k= 10 h/Mpc 



k I h/Mpc 



A^ 

" 1 loop 




-^ 2 ( 9 )G 12 (|k-q|)»^)l} (105) 

where we did not write the time-dependence. At large 
wavenumbers q the integrand Si 2 (q) simplifies to: 

(k.q) 2 

512 (q) ~ , 9 9 [iT > 2lG'ii+i?22Gi2 — R11G21— fll2<J22]-(106) 

Therefore, in the linear regime (or for the steepest-descent 
self-energy E ) where G i} (q; 771, rj 2 ) = e Vl+r,2 P L o(q) the 

two-point correlations Gij factorize and we have: Fig. 18. The density two-point correlation Gu(fc; 771, 772) 

as a function of time 771, for 772 = —1.4 (i.e. z 2 — 3) and 

q >«: Soiia(q) ~ ^L0(g)[i?L;21+i?L;22-i?L;ll-i?L;l 2 ](107) wavenumbers J, = X ( left pancl ) and 1Q x h Mpc -1 ( right 

For a CDM-like power-spectrum P L0 (g) - q" 3 In q at large P and ) ■ We plot the logarithmic power A 2 = 4jt*£Gii. We 
wavenumbers while the linear response RL(q) is constant 
with respect to q. Hence the integration over q yields at 
first sight a logarithmic divergence for q — > 00. In fact, 
we can check from eq. (j2"4"]) that for the linear response the 
terms in the bracket in eq. (|107p actually cancel so that 
the integral is well-defined for the linear two-point func- 
tions (the integrand decreases as 1/q 4 at large q) and it is 
dominated by wavenumbers q ~ k, in agreement with the 
results obtained in sect. 15.11 However, for the 2PI effective 
action method where we must use the non-linear two-point 
functions in eq. (|106p this cancellation no longer holds at 
weakly non-linear scales q^L where they depart from the 
linear predictions. The integral still converges because of 
the decay of the non-linear response at high q but this im- 
plies that for k — > the self-energy term E12 is dominated 
by the contribution of wavenumbers q ~ qNL associated 
with the non-linear transition and not by wavenumbers 
q ~ k. Therefore, E12 does not converge to the steepest- 
descent prediction E 0; i2 at large scales k — > 0. By contrast, 
for other terms the integrand behaves as: 



display the full non- linear power A 2 from eq. (l51|) (solid 
line) and the contribution Af from eq . (|52f (dot-dashed 
line), as well as the usual one-loop perturbative result 
A 2 loop of eq. ([8Tj) (dotted line) and the prediction of the 
steepest-descent method ( "sd" , dotted line) also shown in 
Fig.il 



<?>£;: S , n(q)-S'22(q) 



^RG, S 2 i(q)~^RG. (108) 
q* q A 



The additional power of 1/q ensures that the logarithmic 
divergence of eq. (|107|) does not appear and at large scales 
the self-energy Ey is dominated by the contribution from 
wavenumbers q ~ k hence it matches the results obtained 
within the steepest-descent approach. 

6.2. Correlation G and self-energy II 

Finally, we present in this section our results for the two- 
point correlation G and the self-energy n. We first display 



in Fig. [18] the evolution with time 771 of the density "loga- 
rithmic power" A 2 (fc; 771, 772) defined in eq. ((87j) . As shown 
in the left panel at early times and large scales we recover 
the steepest-descent result since we must also match the 
usual one-loop expansion (|5Tj) . At later times or at smaller 
scales we obtained oscillations which are strongly damped 
at large time separations in contrast with Fig. [8] This 
difference of behaviors of the correlations G between the 
direct steepest-descent method and the 2PI effective ac- 
tion scheme follows from the difference already seen in 
terms of the response R analyzed in sect. 16.11 The right 
panel of Fig.[l8]is particularly interesting as it shows that 
in the highly non-linear regime the two-point correlation 
exhibits a peak at equal times 771 = 772 but whereas both 
the standard one-loop expansion and the steepest-descent 
method yield large oscillations at larger time separations 
the 2PI effective action approach leads to a strong damp- 
ing. Thus, this is the only method (among those presented 
in this paper) which yields (at one-loop order) a damping 
of correlations at small scales and large time separations 
which is expected to reflect the qualitative behavior of the 
exact gravitational dynamics. 

Next, we show in Figs. [19l [20] the logarithmic power 
A 2 (fc; 771, 772) at equal redshifts z\ — z 2 = and Z\ = z 2 = 
3 as a function of wavenumber k. We find again that our 
results match the steepest-descent prediction as well as the 
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Fig. 19. The logarithmic power A 2 (fc) at redshift z = 0, 
that is at equal times z\ = z 2 — 0. We display the full non- 
linear power A 2 from eq. (|51[) (solid line) and the contribu- 
tion A 2 from eg . ([52f (dot-dashed line), as well as the linear 
power A 2 (dashed line), the usual one- loop perturbative 
result A 2 loop of eq. (|8Tj) (dotted line), the prediction of the 
steepest-descent method ( "sd" , dotted line) also shown in 
Fig. [5] and the fit A 2 b from numerical simulations (Smith 
et al. 2003). All quantities (including A 2 ) are positive. 
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Fig. 20. The logarithmic power A 2 (fc) at redshift z = 3, 
that is at equal times Z\ = Z2 = 3. 

usual one-loop power (f5Tj) at large scales. At small scales, 
despite the damping at unequal times shown in Fig. [TJ] 
for the response function and in Fig. [TH] for the corre- 
lation function we obtain a steady growth of the power 
A 2 (fc) in between the linear prediction A 2 L and the usual 
one- loop prediction A 2 loop . Besides, the agreement with 
the results from numerical simulations is better than for 
the stcepcst-dcscent prediction. Of course, because of the 



Fig. 21. The logarithmic power A 2 (k) at unequal times 
(m = 0,772 = -1.4) (i.e. ^=0,22=3). 

damping of the response R at large time separations an- 
alyzed above we can check that the contribution A 2 from 
eq. (|52|) associated with the transport of the initial fluctua- 
tions becomes negligible in the non-linear regime (contrary 
to what happens within the steepest-descent approach, see 
Figs. O [TO]). However, the continuous generation of fluctu- 
ations, described by the self-energy II in eq. (|5Tj) and asso- 
ciated with the non-linearity of the dynamics, sustains a 
high level of density fluctuations A 2 (fc) into the non-linear 
regime. The latter seen at equal times 7/1 = 772 = f) are re- 
lated to the values of n(fc; r]' 2 ) at nearly equal times 
V'l ^V2-V- 

We show in Fig.[2T]the logarithmic power A 2 (fc; 771, 772) 
at unequal times as a function of wavenumber k. We find 
again that our result matches the steepest-descent pre- 
diction as well as the usual one-loop power (|5Tj) at large 
scales. However, the correlation now quickly decays at 
small scales in the non-linear regime. This follows from 
the damping at large time separations analyzed above for 
the response function which leads to a decorrelation at 
small scales and large time separations. 

We display in Fig. [55] the self-energy II as a function 
of forward time 771. We do not show the results of the 
steepest-descent method which follow a mere exponential 
growth from eq. ([65|) At weakly non-linear scales (k ~ 1/t, 
Mpc -1 ) we recover at early times the exponential growth 
of the steepest-descent prediction whereas at late times we 
obtain a fast decay. At highly non-linear scales (k = lOh 
Mpc -1 ) we obtain a peak at equal times 771 = 7/2 and 
strongly damped oscillations at large time separations. 
This is clearly consistent with the results obtained for the 
correlation Gn displayed in Fig. [T5I 

Next, we show in Fig. [23] the self-energy n at equal 
redshifts Z\ = zi = as a function of wavenumber 
k. We again recover the steepest-descent prediction at 
large scales. However, at small scales we now obtain a 




Fig. 22. The self-energy ±l(fc; 771,772) as a function of for- 
ward time 771, for 772 = —1.4 (i.e. z 2 = 3) and wavenumbers 
k = 1 and 10 x h Mpc" 1 . 



1000 











n 


/ /-_ 




Z 1 = Z 2=° / 






y — 




/ <,?/ 
/,?/ 


sd ^ 




/>'/ / 
// /' 

// 1 

f" ' 

i 1 1 1/. 1 





0.01 0.1 1 10 



k [h Mpc- 1 ] 

Fig. 23. The self-energy H(k) as a function of wavenum- 
ber k at equal redshifts z\ = z 2 = 0. We also plot the 
prediction of the steepest-descent method (curves labeled 
"sd"). 



steady growth whereas the steepest-descent method yields 
a smooth decline. Therefore, in agreement with Figs. I19[ 
l20l we find that the coupled system of equations (|46l) 
(with Go replaced by G) and (f5Tj) entails at equal times 
a stronger growth of the correlation G and of the self- 
energy II into the non-linear regime as compared with 
the non-coupled equations associated with the steepest- 
descent method. 

Finally, we display in Fig. [24] the self-energy II at un- 
equal redshifts z± = 0, z 2 = 3 as a function of wavenumber 
k. In agreement with our results for other two-point func- 



Fig. 24. The self-energy II(fc) as a function of wavenumber 
k at unequal redshifts z\ = 0, Z2 = 3. We also plot the 
prediction of the steepest-descent method (curves labeled 
"sd"). 

tions we recover the steepest-descent prediction at large 
scales and decaying oscillations at small scales. 

7. Simplified response approximations 

The results obtained in the previous sections suggest that 
a convenient approximation would be to use a simple an- 
alytical form for the response function R and to compute 
the correlation G and the self-energy II from eqs. (|46|) , (f5l"j) . 
This avoids the computation of the self-energy term S 
and the numerical integration of the integro-differential 
eq. (|43p which requires rather small time-steps. By con- 
trast, eq. ([5Tj) is a Volterra integral equation of the sec- 
ond kind which is usually better behaved. We shall in- 
vestigate two such approximations associated with the di- 
rect steepest-descent approach and the 2PI effective action 
scheme. 

7.1. Simplified response for the steepest-descent 
scheme 

Within the framework of the direct steepest-descent 
scheme studied in sect. [5] we found that the response func- 
tion exhibits in the non-linear regime a series of oscilla- 
tions with an envelope given by the linear theory predic- 
tion. Therefore, following eq. (|80p we consider the approx- 
imate response i? app defined as: 

R app (x ll x 2 ) = R L {xi 1 x 2 )cos[uj{k)(ai - a 2 )], (109) 

where Rl{= Ro) is the linear response given in ea. ([M|) . 
a = e'' is the scale factor and u>{k) is given by eg. ((77)) . 

Here it is interesting to see how expression (|109|) may 
be interpreted. In the linear regime the density and ve- 
locity fluctuations ipL are merely amplified by a time- 
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dependent factor (proportional to the scale-factor) which 
implies in particular that the physics is exactly local: 
ipL (x, 77) only depends on the initial conditions (and a pos- 
sible external noise) at the same location x. Thus in real 
space we have -Rl(x 1 ,x 2 ) cx 6jj (xi — x 2 ) and in Fourier 
space, using the factorization ([50)1 associated with statis- 
tical homogeneity, -Ri(fc) is constant with respect to k, 
in agreement with cq. (l24]) . On the other hand, the cosine 
dependence of eq. (|109p yields in real space: 



Ra,p P (xi,x 2 ) = RL(ai,a 2 )V(x-i - x 2 ; a%, a 2 ) 



(110) 



where i?i(ai,a 2 ) is the time-dependent factor of the lin- 
ear response (|24l) and we introduced the distribution 
■p(r;ai,a 2 ) obtained from eq. (|109[) as: 

7>(r; oi, a 2 ) = — —S' D (r - d) with d = (ai - a 2 )r (111) 
Airr 



where r — Irl and: 



To 



47T f 

— I dqP L0 (q) 



1/2 



(112) 



Thus, eqs. (|109|) - (|112|) may be interpreted as the displace- 
ment of the linear fluctuations over a random distance r 
with a probability distribution V{r) (whereas the linear 
regime obviously corresponds to the distribution (5o(r)). 
This non-linear operation corresponds to some diffusion 
of the initial conditions (over angles of the shift r) but 
the singular distribution over the length |r| leads to the 
mere cosine tail in Fourier space instead of a strong decay. 
This can be contrasted to the Zeldovich approximation 
(Zeldovich 1970) where the trajectory of particles with ini- 
tial comoving location q is written as x(q, a) = q + r with 
V q .r = — <5i(q, a). Thus the displacement field is Gaussian 
which leads to a Gaussian decay at high k (Crocce & 
Scoccimarro 2006a). However, we must emphasize that we 
only write eq. (|110p for comparison but this is not neces- 
sarily the best way to understand the result of the direct- 
steepest approach (at this stage expression (|110[) is only 
a mere re- writing of eq. (|109p k In particular, there is no 
reason why the non-linear response should be factorizcd 
as in eq. ()110p as Rl times a displacement term. 

Next, we use the steepest-descent prescription (|65j) - 
(f66|) for the self-energy II and we compute numerically the 
correlation G from the explicit expression (|5ip . In fact, 
using eq. (|65p the integrals over time may be performed 
analytically which yields: 



G app (fc;?7i,f?2) 



with: 



Gh{k] 771 , 7? 2 ) cos(wai) cos(wa 2 ) 
+a 2 l a\R(<jjai).IiQ(k).R(ua 2 ) T 



(113) 



R{ua) = f+(Lua)R+ + /_(wa)i^. 

The matrices R^ are given in eqs. 
the functions: 



„ . . f 1 , , ,„ sin(a;) 
f+(x) = / dt cos[x(l - t)] = — — 



(114) 

and we introduced 
(115) 
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Fig. 25. The logarithmic power A 2 (fc) at redshift z = 0. 
We display the ratio A 2 (fc)/A| / (fc), where A 2 L is the linear 
power, for the usual one-loop perturbative result A 2 loop 
(dotted line), a fit A 2 b (upper solid line) from numerical 
simulations (Smith et al. 2003), the result of the steepest- 
descent method ( "sd" , dotted line) shown in Fig. [3] and 
the approximate power A 2 pp (solid line) obtained from 
the approximation (|109p . 
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Fig. 26. The ratio A 2 (fc)/A|(fc) as in Fig 
shift 2 = 3. 



lbut at red- 



and: 



f-(x) = / dtt 5/2 cos[x(l-t)} 
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cos(x)S{yj2x/*K) - sm.{x)C{^/2xJ^) , (116) 



where C(x) and S(x) are the Fresnel integrals. At early 
times or large scales the first term dominates in eq. (IH3|l 
and we recover the linear regime G — > Gl- At late times 
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and small scales we must take into account both terms 
and for large uja,i we obtain: 

Luai 3> 1 : G(fc; 0,1,0,2) ~ a\a,2^Gw{k) cos(wai) cos(wa2) 

+R+.n (k).Rr sin( " ai y n( " fl2) }. (117) 

Thus, we can see that at late times the non-linear corre- 
lation grows as G(fc; 01,02) ~ 0102 with oscillations for 
0,1 7^ &2- This agrees with the results presented in sec- 
tion [Ol 

We show our results for the equal-time density correla- 
tion G11 (i-e. the power A 2 (A:)) in Figs. I251 !2fJl at redshifts 
z = and z = 3. We display the ratio to the linear power 
A 2 (fc)/A 2 (fc) to zoom over the weakly non- linear region 
and we also plot the exact one-loop steepest-descent pre- 
diction ("sd") shown in Fig. [5] for comparison. Figs. 1251 
[251 clearly show that the exact one-loop steepest-descent 
prediction matches the standard one-loop perturbative re- 
sult at large scales. Moreover, Fig. [25] suggests that at 
z = the fit from numerical simulations (Smith et al. 
2003) may be off by 10% at k ~ 0.3/i Mpc" 1 . The power 
A 2 pp follows roughly the behavior of the "exact" one-loop 
steepest-descent prediction but it is smaller by up to 15% 
and does not match very well the one-loop perturbative 
result. This failure at intermediate scales is not really sur- 
prising since eq. ()80|) was derived in the large-fc limit using 
the asymptotic behaviors (163l) - (l64]l . However, as can be 
seen in Fig. Q]they are only reached for k > lOh Mpc -1 . 
Therefore, the approximation (]109[) is not sufficient to ob- 
tain accurate results over weakly non-linear scales. This 
shows that the results obtained for the two-point corre- 
lation are rather sensitive to the approximations used for 
the response R. In fact, since the direct steepest-descent 
method is quite simple and easy to compute the approx- 
imation (| 1 09|) is not very useful for precise quantitative 
predictions and it is best to use the rigorous method pre- 
sented in sect. [5] 

7.2. Simplified response for the 2PI effective action 
scheme 

Within the 2PI effective action method we found in sect. [6] 
that the response function exhibits oscillations as for the 
steepest-descent prediction but their amplitude decays at 
large times. Following the analysis of the simple model 
(|95p and eq. (1104p we consider the approximate response 
-Rapp defined as: 



-Rap P (xi, x 2 ) = Rl(xi,x 2 



^ Ji[2u(k)(fli - 02)] 
w(fc)(oi - a 2 ) 



(118) 



Of course, this expression can again be written in the form 
(fTT0|) in real space, with a distribution V(r) which is now 
given by: 



where the top-hat factor 0(r < 2d) obeys 6 = 1 for r < 2d 
and 9 = for r > 2d. The decay of the response func- 
tion at high k corresponds to the smoother distribution 
()119l) as compared with eq. (Illlj) , but the discontinuity at 
r = 2d leads to a slow power-law damping instead of an 
exponential decay. Again, ea. (|119p is a mere rewriting of 
ea. (|118p and we do not claim here that the distribution 
Vir) is anything more than a phenomenological tool. 



Then, we use eq. (|46)l (with Go replaced by G) and 
eq. (f5"Tj) to compute the self-energy II and the correlation 
G, following the procedure described in sect. [6] 
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Fig. 27. The logarithmic power A 2 (k) at redshift z = Q. 
We display the ratio A 2 (fc)/A| / (fc), where A 2 is the linear 
power, for the usual one-loop perturbative result A 2 loop 
(dotted line), a fit A 2 b (upper solid line) from numerical 
simulations (Smith et al. 2003), the result of the 2PI effec- 
tive action method ( "2PI" , dotted line) shown in Fig. [T5] 
and the approximate power A 2 pp (solid line) obtained 
from the approximation (|118p . 



We show our results for the power A 2 (fc)) in Figs. [27l 
l28l at redshifts z = and z = 3. We again display the 
ratio to the linear power A 2 (/c)/A 2 (fc) and we also plot 
the exact 2PI effective action prediction ("2PI") shown 
in Fig. [TH] for comparison. Figs. [27] [25] again clearly show 
that the exact 2PI effective action prediction matches the 
standard one-loop perturbative result at large scales. In 
a fashion similar to the approximation (j!09[) we find that 
the power A 2 pp follows roughly the behavior of the exact 
one-loop 2PI effective action prediction but it is smaller 
by up to 15% and does not match very well the one-loop 
perturbative result. This shows again that the two-point 
correlation is rather sensitive to the value of the response 
function. 
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Fig. 28. The ratio A 2 (fc)/A|(fc) as in Fig. W\ but at red- 
shift z — 3. 

7.3. Mixed approaches 

Finally, we consider in this section mixed approaches 
which use elements from both the steepest-descent and 
2PI effective action schemes. The idea is to simplify the 
computation involved in the 2PI effective action method 
by separating the computations of the pairs (S, R) and 
(H, G). That is, we first compute the self-energy S and 
the response R from the coupled eqs. (f4"3"]) . (|45[) . where we 
use Go and R in the expression ([45]) for £ (the steepest- 
descent scheme uses Go and Rq whereas the 2PI effective 
action scheme uses G and R). Thus, the pair (£, R) be- 
comes independent of (II, G) and can be computed in a 
first step. Since we keep the explicit dependence on R in 
the self-energy £ the evolution equation (|43f for the re- 
sponse function is no longer linear but quadratic, hence we 
expect to recover the damping discussed in sect. 16.11 (see 
the analysis of eq. ([95|) ). We display our results in Fig. [29] 
which shows that the non-linear response exhibits indeed 
a damping close to the one obtained in Fig. [TJ] for the full 
2PI effective action scheme. 

Next, we compute the pair (II, G) from the equations 
(|4"6"|) and ([51]) . Here we can consider two procedures as 
we may use either Go or G in the computation of II. 
Clearly, these mixed schemes still agree with the usual one- 
loop perturbative results. These two-step approaches are 
slightly faster than the full 2PI effective action method for 
numerical computations. Indeed, as explained in sect. [51 
within the 2PI effective action approach the equations 
obtained for two-point functions and self-energies form a 
coupled system of non-linear equations. This leads to an 
iterative scheme for their numerical computation (which 
converges within ~ 7 steps at worst for the scales we con- 
sider here). On the other hand, at each step the computa- 
tion of the correlation G takes most of the computer time 
(especially at late times) because of the double integral 
over times in eq. ([5T|) but it converges faster (at fixed R) 




k [h Mpc- 1 ] 



Fig. 29. The non-linear response function R(k; rji, 772) as a 
function of wavenumber k, at times 771 = 0,772 = —1-4, for 
the mixed approaches where Go is used in the computation 
of S. 

than the response R (at fixed G). Therefore, it saves time 
to separate the computations of R (many fast iteration 
steps) and G (few long iteration steps). 
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Fig. 30. The logarithmic power A 2 (fc) at redshift z = Q, 
that is at equal times z\ — Z2 = 0. We display the lin- 
ear power A 2 (dashed line), the usual one- loop perturba- 
tive result A 2 loop (dotted line) , a fit A 2 b (lower solid line) 
from numerical simulations (Smith et al. 2003), the mixed 
approach A 2 ^ based on IIo (dashed line) and the mixed 
approach A^ where II is coupled to G (upper solid line). 

We display in Figs. [3011331 our results for the power 
A^fc) at redshifts z = and z = 3 for both schemes. The 
dashed curve Afj shows the power obtained when the 
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Fig. 31. The logarithmic power A 2 (fc) at redshift z = 3, 
that is at equal times z\ = z-i = 3. 
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Fig. 32. The ratio A 2 (fc)/A 2 (fc) at redshift z = of the 
logarithmic power A 2 to the linear power A 2 . We display 
the ratio A 2 /A| for the usual one-loop perturbative re- 
sult A 2 loop (dotted line), a fit A 2 b (lower solid line) from 
numerical simulations (Smith et al. 2003), the mixed ap- 
proach A|[ o based on IIo (dashed line) and the mixed ap- 
proach A^ where II is coupled to G (upper solid line). 

self-energy term LTo is used (as in the steepest-descent ap- 
proach) so that the correlation G is directly given by the 
explicit expression (|5ip . The solid curve A 2 ^ corresponds 
to the case where we use the non-linear correlation G in 
the expression ([46]) of the self-energy, so that eqs. (|46l) 
and (|51[) are a coupled non-linear system which we solve 
iteratively (as in the 2PI effective action scheme) . We can 
see that both methods exhibit a similar behavior until 
A 2 (k)/A 2 L (k) ~ 1.4. Farther into the non-linear regime 
the "LTo-scheme" breaks from the non-linear growth of 
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Fig. 33. The ratio A 2 (fc)/A|(fc) as in Fig. [32] but at red- 
shift z — 3. 

the power-spectrum and shows oscillations which con- 
verge close to the linear power, in a manner reminiscent 
of the steepest-descent approximation, whereas the cou- 
pled "Il-scheme" follows this non-linear growth. On the 
other hand, in the weakly non-linear regime displayed in 
Figs. I32ll33l the predicted power A 2 (k) is suppressed as 
compared with the usual one-loop result. Depending on 
the redshift this can improve or worsen the global agree- 
ment with numerical simulations. However, Fig. [33] sug- 
gests that when the usual one-loop result shows a better 
global agreement with simulations it is a mere coincidence 
as the simulation result separates last from the large- 
N predictions (and the closer agreement with the usual 
one-loop result at smaller scales has no theoretical basis). 
Nevertheless, the limited accuracy of the numerical simu- 
lations prevents one from drawing definite conclusions. We 
can note in Fig. [30] that in the highly non- linear regime 
the approximation A^ grows larger than for the full 2PI 
effective action prediction displayed in Fig. [T9] Indeed, in 
the latter case the non-linear response is further damped 
by the strong growth of the two-point correlation G (this 
would correspond to a larger a whence u> in the simple 
model (|104jl ) and this "negative feedback" within the cou- 
pled system R, G leads to a smaller G as compared with 
the approximation A^ where this feedback is neglected. 
In fact, at very high k (beyond the range shown in Fig. 130]) 
the power A^ appears to diverge. It is not clear whether 
this is due to the finite resolution of the numerical com- 
putation or to a true divergence. However, for practical 
purposes this is not a serious shortcoming since it occurs 
at small scales beyond the range of validity of this ap- 
proach (one would need to go beyond one-loop order and 
probably beyond the hydrodynamical approach as shell- 
crossing can no longer be neglected). 

The good agreement in the weakly non-linear regime of 
these various schemes enables one to evaluate their range 
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of validity without the need to compare their prediction 
with N— body simulations, which could be of great prac- 
tical interest. In fact, Fig. I3"2l suggests that this procedure 
can provide a better accuracy than numerical simulations 
in this regime, where the latter may still be off by 10%. 

8. Application to a ACDM universe 

We finally describe in this section how the results ob- 
tained in the previous sections can be applied to more 
general cosmologies, such as a ACDM universe. As seen 
in sect. !2.1l our formalism applies equally well to any cos- 
mology. We only need to use the relevant matrix 0, that 
is to take into account the time dependence of the ratio 
^m// 2 in ea. lfTB")) . and to recall that r\ is the logarithm of 
the linear growth factor, see eq.©. In practice, a widely 
used approximation is to take Q m / f 2 — 1 so that the 
results obtained for the critical-density universe directly 
apply to the case A ^ and the only change is the time- 
redshift relation r\ <-> z. We shall use this simple approx- 
imation here and compare our results to numerical sim- 
ulations. Thus, we consider the cosmological parameters 
n m = 0.3, 0\ = 0.7, er 8 = 0.9 and H = 70 km.s-kMpc" 1 
as in Smith et al.(2003). 
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Fig. 34. The non-linear response function 771 , 772) 
given by the direct steepest-descent method as a func- 
tion of wavenumber k for a ACDM universe, at redshifts 
Z\ = 0, Z2 = 3. The horizontal lines are the linear response 
Rl of eq.(|24[) which is independent of k. 

We first plot in Figs. I34|35l the non-linear response 
functions we obtain within the one-loop steepest-descent 
and 2PI effective action methods. We can check that we 
recover the behavior obtained for the critical-density uni- 
verse displayed in Figs. 171141 The steepest-descent ap- 
proach yields again oscillations in the non-linear regime 
with an amplitude given by the linear response at high 
k whereas the 2PI effective action method gives damped 
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Fig. 35. The non-linear response function R(k; 771 , 7/2) 
given by the 2PI effective action method as a function 
of wavenumber k for a ACDM universe, at redshifts z\ = 
0, Z2 — 3. We also plot the first half-oscillation of the re- 
sponse obtained from the direct steepest-descent method 
(curves labeled "sd" ) which was shown in Fig. l34l 



oscillations. The transition to the non-linear regime oc- 
curs at slightly smaller wavenumber because of the larger 
normalization as of the linear power-spectrum. 
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Fig. 36. The logarithmic power A 2 (k) at redshift z = 
for a ACDM universe. We display the linear power 
A 2 (dashed line), the usual one- loop perturbative result 
A 2 loop (dotted line), a fit A 2 b (lower solid line) from 
numerical simulations (Smith et al. 2003), the steepest- 
descent prediction (dotted line "sd" ) and the 2PI effective 
action result (solid line "2PI" ) . 
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Fig. 37. The logarithmic power A 2 (k) at redshift z = 3 Fig. 39. The ratio A 2 (k)/A 2 L (k) as in Fig.[38]but at red- 
fbr a ACDM universe. shift z = 3. 




0.05 0.1 0.5 

k [h Mpc- 1 ] 



Fig. 38. The ratio A 2 (fc)/A|(fc) at redshift z = of the 
logarithmic power A 2 to the linear power A\ for a ACDM 
universe. We display the ratio A 2 /A| for the usual one- 
loop perturbative result A 2 loop (dotted line), a fit A 2 b 
(lower solid line) from numerical simulations (Smith et al. 
2003), the steepest-descent prediction (dashed line "sd") 
and the 2PI effective action result (solid line "2PI" ) . 

We show in Figs. I3l)ll3l?l our results for the two-point 
correlation Gu, that is the non-linear density power spec- 
trum. We again recover the behavior analyzed in previous 
sections for the critical-density universe. Both large- N ex- 
pansion schemes agree with the standard one-loop expan- 
sion at large scales. At small scales the steepest-descent 
result goes back to values close to the linear power A 2 
whereas the 2PI effective action result keeps growing in a 
manner similar to the usual perturbative result and nu- 
merical simulations. The degree of agreement between the 



N— body simulations and either one of the standard and 
2PI results depends on redshift, that is on the slope of the 
linear power-spectrum at the relevant scales. 

Thus, the large- N expansions developed in this article 
apply in the same manner for different cosmologies and 
exhibit the same properties. 

9. Conclusion 

In this article we have studied the predictions at one-loop 
order of large- N expansions for the two-point correlations 
associated with the formation of large-scale structures in 
the expanding universe. Focussing on weakly non-linear 
scales we have used the hydrodynamical equations of mo- 
tion as a starting point. Then, we have recalled the path- 
integral formalism which describes the statistical proper- 
ties of the system, assuming Gaussian initial conditions, 
and we have presented two possible large-./V expansions. 
Next, we have described in details the numerical predic- 
tions obtained for a SCDM cosmology. 

First, we have presented a direct steepest-descent ap- 
proach where the self-energy is expressed in terms of the 
linear response and correlation. This allows simple com- 
putations as the time-dependence of the self-energy can 
be factorized (for a critical density universe). Moreover, 
the non-linear response R can be computed directly from 
this self-energy as it is not coupled to the non-linear cor- 
relation G and it obeys a simple linear equation. This also 
allows a detailed analysis of its behavior. We have found 
that at one-loop order this approach yields a non-linear 
response which exhibits strong oscillations in the non- 
linear regime (small scales or late times) with an ampli- 
tude which is given by the linear response. Therefore, the 
response does not decay but only exhibits increasingly fast 
oscillations of the form i?£COs[(ai — d2)k/k*\. However, 
once the response is integrated over time such oscillations 
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lead to some damping as compared with the linear re- 
sponse. Then, we have computed the non-linear correla- 
tion which can be explicitly written in terms of this non- 
linear response which transports forward over time the 
initial density and velocity fluctuations as well as the fluc- 
tuations generated at all previous times by non-linear cou- 
plings. We have checked that in the quasi- linear regime the 
non-linear correlation agrees with the usual one-loop pre- 
diction (obtained from the standard perturbativc expan- 
sion over powers of the initial fluctuations) . In the highly 
non-linear regime the correlation separates from the steep 
growth displayed by the standard one-loop result and con- 
verges back to the linear amplitude. Unfortunately, this 
behavior does not improve the agreement with N— body 
numerical simulations. Nevertheless, it suggests that this 
large- N expansion is well-behaved as the two-point func- 
tions do not "explode" in the highly non-linear regime. 
Contrary to the standard perturbative expansion, where 
higher-order terms grow as increasingly large powers of 
k and a(t), the higher-order terms obtained within this 
method may remain of finite amplitude. 

Secondly, we have described a 2PI effective action ap- 
proach where all two-point functions are coupled. The nu- 
merical computation requires an iterative procedure but 
thanks to causality the equations can be solved by inte- 
grating forward over time with the help of a few iterations 
at each time-step (keeping an implicit scheme in the non- 
linear intcgro-diffcrcntial equations). Then, we have found 
that the one-loop non-linear response again exhibits fast 
oscillations in the non-linear regime but its amplitude now 
shows a fast decay. We have explained how this damping 
is produced by the non-linearity of the evolution equation 
obeyed by the response function. This yields a power-law 
decay such as i?z,Ji[2(ai — a,2)k/k*\/[(ai — a,2)k/k*\. In the 
quasi-linear regime the non-linear correlation agrees again 
with the usual one-loop prediction but in the highly non- 
linear regime it no longer converges back to the linear 
power. It rather shows a moderate growth intermediate 
between the linear prediction and the result of numerical 
simulations, with an amplitude which depends on the red- 
shift (whence on the slope of the linear power-spectrum). 
Therefore, this method appears to be superior to the di- 
rect steepest-descent approach which could not recover 
any non-linear damping (at one-loop order) but it requires 
more complex numerical computations. 

Thirdly, we have described other approximation 
schemes built from these two approaches. Thus, we have 
found that using a simple explicit approximation for the 
non-linear response within cither large-A scheme gives 
correct qualitative results but does not manage to repro- 
duce accurately the quantitative results obtained from the 
full large-A approaches. Next, we have investigated the 
approximations obtained by separating the computations 
of the response and the correlation (but keeping a non- 
linear dynamics for the response). As expected we have 
found that this yields a non-linear response which closely 
follows the damping obtained by the full 2PI effective ac- 



tion approach and it gives good results (as compared with 
the full large- N schemes) for the correlation. 

Finally, we have shown that these results also apply to 
more general cosmologies such as a ACDM universe, since 
our formalism extends to any expansion history of the 
Universe provided we use the correct time-redshift relation 
and linear growth factors. 

Thus, we have shown in this article that large- A ex- 
pansions provide an interesting systematic approach to the 
formation of large-scale structures in the expanding uni- 
verse. They already show at one-loop order some damp- 
ing in the non-linear regime for the response function and 
the amplitude of the two-point correlation remains well- 
behaved. Unfortunately, for practical purposes the accu- 
racy provided by these large-A expansions at one-loop or- 
der is not significantly better than the usual perturbative 
result for weakly non-linear scales. Nevertheless, the im- 
provement may become greater at higher orders since the 
expansion is likely to be better behaved, but this requires 
rather complex computations which are beyond the scope 
of this paper. 

On the other hand, one of the goals of this article 
was to describe a different approach to the problem of 
non-linear gravitational clustering than the standard ex- 
pansion schemes, namely a path-integral formalism. An 
important feature of this approach is that the statistical 
nature of the problem is already included in the defini- 
tion of the system, that is the action S[ip, A] describes 
both the equations of motion (here in the hydrodynami- 
cal limit) and the average over Gaussian initial conditions. 
This also means that one directly works with correlation 
functions (such as the power spectrum). By contrast, in 
usual expansion schemes (or in A-body numerical simu- 
lations) one first works with the density field associated 
with a specific initial condition, that is one tries to ex- 
press the non-linear density field in terms of the linear 
density field (e.g. as a truncated power series) and next 
performs the average over the initial condition. In princi- 
ples, this formalism can present several advantages. First, 
the statistical quantities such as the power-spectrum are 
precisely the objects of interest for practical purposes. 
Secondly, they exhibit symmetries (such as homogeneity 
and isotropy) which are not satisfied by peculiar realiza- 
tions of the initial conditions. Moreover, they are smooth 
functions (such as power laws) rather than singular distri- 
butions. This could facilitate the numerical computations 
and the building of various approximations. 

Of course, one can recover the standard perturbative 
results by expanding over the cubic part of the action, 
as recalled in this paper and discussed in more details in 
Valageas (2004). However, the advantage of this formalism 
is that it may serve as a starting point to other approx- 
imation schemes by applying the methods of statistical 
mechanics and field theory. Thus, we have described in 
details in this paper two large-A expansion schemes. In 
addition to the interesting properties of these two meth- 
ods discussed above, we can note that it is useful to 
have several rigorous expansion schemes which only differ 
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by higher-order terms than the truncation order. Indeed, 
since one can expect that they should be correct up to the 
scale where they start to depart from one another, this 
should allow one to estimate their range of validity with- 
out computing explicitly higher-order terms or perform- 
ing A^-body simulations. In fact, in the weakly non- linear 
regime where all expansion schemes agree this procedure 
should provide a better accuracy than iV-body simulations 
which may still be inaccurate by up to 10%. Farther into 
the non-linear regime, it seems that the usual one-loop re- 
sult and the 2PI effective action scheme give better results 
than the direct steepest-descent method. On the other 
hand, one can hope that other methods than these two 
expansion schemes could be applied to the action S[ip, A]. 

In regard to recent works we can note that 
the Schwinger-Dyson equations obtained in Crocce & 
Scoccimarro (2006a) can also be obtained from the ap- 
proach described in this paper, provided we use the equa- 
tions of motion in their integral form (that is we first in- 
tegrate the time-derivative) , as also described in Valageas 
(2001). Whereas Crocce & Scoccimarro (2006a) used a 
diagrammatic technique and the resummation of all di- 
agrams gives back the usual Schwinger-Dyson equations 
the path-integral formalism directly gives these equations 
by expanding the action (or its Legendre transform) about 
a saddle-point. On the other hand, the diagrammatic ap- 
proach also provides a direct expression of the response 
function R(r)i,T)2) (in the limit 772 — > —00) as a series of 
diagrams (which could be recovered here by writing the 
self-energy £ as a series of diagrams). Thus, Crocce & 
Scoccimarro (2006a) managed to resum these diagrams in 
the high-fc limit and to obtain the asymptotic behavior 
of i?(?7i,?72 — * —00). It is not clear yet how this result 
could be directly obtained from the path-integral formal- 
ism without expanding over diagrams. 

Among the points which deserve further studies, it 
would be interesting to apply these approaches to higher- 
order correlations beyond the two-point functions inves- 
tigated here. On the other hand, the analysis presented 
in this article is based on the hydrodynamical equations 
of motion which break down at small scales beyond shell- 
crossing. It is possible to apply these large- N expansions 
to the Vlasov equation (Valageas 2004) but the numer- 
ical computations would be significantly more complex 
(since one needs to add velocities to space coordinates). 
Nevertheless, the results obtained within the hydrody- 
namical framework (damping and well-behaved non-linear 
asymptotics) can be expected to remain valid in the col- 
lisionless case and studies such as the present one may 
serve as a first step towards an application to the Vlasov 
dynamics. Such large- N expansions could also be applied 
to simpler effective dynamics which attempt to go be- 
yond shell-crossing (based for instance on a Schroedinger 
equation, Widrow & Kaiser 1993). These works would be 
of great practical interest as several cosmological probes 
(such as weak lensing surveys and measures of the baryon 
acoustic oscillations) which aim to measure the recent ex- 
pansion history of the Universe in order to constrain the 



dark energy equation of state (as well as other cosmologi- 
cal parameters) are very sensitive to the weakly non-linear 
scales where non-linear effects can no longer be neglected. 
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